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Abstract 


The  thrust  of  this  research  program  has  been  the  improvement  of  our  capa¬ 
bilities  for  analyzing  stability  and  transition  in  compressible  boundary-layers 
as  they  appear  in  technologically  important  flows.  Examples  of  such  flows  are 
swept  wings  of  commercial  airplanes  or  the  highly  cambered  blades  of  gas 
turbines.  We  have  extended  the  parabolized  stability  equations  (PSE)  for 
these  situations  and  developed  methods  for  solving  these  equations  in  distur¬ 
bance  environments  reaching  from  the  low  atmospheric  levels  to  the  volatile 
high  levels  in  gas  turbines.  The  extension  required  both  a  basic  study  on  the 
stability  analysis  of  3D  flows  and  formulating  the  equations  in  general  curvi¬ 
linear  coordinates.  Proper  choices  have  been  made  to  minimize  the  effect 
of  the  parabolization.  Major  efforts  have  been  spent  on  characterizing  the 
disturbance  environment.  More  recently,  emphasis  has  been  directed  toward 
the  control  of  stability  and  transition  using  neural  networks. 
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1  Objectives 

This  research  program  aimed  at  developing  and  applying  theoretical,  numer¬ 
ical,  and  graphical  tools  for  the  quantitative  description  and  deeper  under¬ 
standing  of  stability  and  transition  in  technologically  important  flows.  Much 
of  the  effort  has  been  directed  toward  engineering  methods  for  advanced 
design.  Toward  these  goals,  we  have  worked  in  the  following  areas: 

(1)  Optimum  formulations  and  algorithms  for  solving  the  PSE  in  general 
curvilinear  coordinates. 

(2)  High-order  Hermitian  finite-difference  methods. 

(3)  Stability  theory  for  3D  boundary  layers. 

(4)  Basic-flow  calculations  with  Euler  and  boundary-layer  codes.  Accuracy 
requirements. 

(5)  Transition  mechanisms. 

(6)  Characterization  of  the  disturbance  environment.  Receptivity  studies. 

(7)  Transient  growth  of  disturbances  in  stable  regions. 

(8)  Modeling  of  nonlinear  phenomena  with  neural  networks.  Aspects  of 
flow  control. 

2  Achievements 

The  work  in  all  areas  has  considerably  progressed.  We  have  developed  two 
versatile  PSE  codes.  The  first  code  is  very  flexible  to  adapt  to  2D  or  3D  con- 
vectively  unstable  flows.  The  code  is  restricted  to  simple  (typically  Caxtesian) 
coordinate  systems  and  primarily  serves  for  producing  benchmark  results  and 
exploring  new  concepts  using  exact  or  similarity  solutions  for  the  basic  flow. 
The  second  code  has  been  developed  in  cooperation  with  DynaFlow,  Inc.  and 
is  written  specifically  for  compressible  boundary  layers  in  general  curvilinear 
coordinates  to  enable  studies  of  flows  in  realistic  geometries.  This  code  is  the 
basis  for  specific  applications  e.g.  to  transition  and  heat  transfer  analysis  in 
gas  turbines.  The  work  with  this  code  is  supported  by  various  other  codes 


3 


to  compute  the  inviscid  flow  over  wings  or  through  cascades  and  to  obtain 
the  boundary-layer  flow  from  the  inviscid  surface-pressure  distribution.  The 
work  under  this  contract  has  involved  various  graduate  students  one  of  whom 
will  complete  his  Ph.D.  in  the  near  future.  The  difficulty  of  the  research  ar¬ 
eas  listed  above  has  required,  however,  a  major  involvement  of  the  principal 
investigator.  Some  details  on  progress  in  the  various  areas  of  interest  are 
reported  below. 

2.1  PSE  in  General  Curvilinear  Coordinates 

Previous  work  has  either  neglected  or  only  partially  accounted  for  the  cur¬ 
vature  effects  on  basic  flow  and  stability  characteristics.  While  the  step-by- 
step  inclusion  of  curvature  in  previous  work  has  provided  the  opportunity 
for  follow-up  publications  of  minor  value,  we  have  decided  to  account  for  the 
curvature  terms  completely  and  from  the  beginning.  In  a  first  approach,  we 
have  derived  the  nonlinear  incompressible  stability  equations  including  trans¬ 
verse  and  longitudinal  curvature  as  functions  of  the  distance  from  the  wall. 
The  linearized  equations  have  been  coded  for  use  with  both  spectral  method 
or  compact  finite-difference  method.  The  coded  insert  files  are  consistent 
with  the  file  format  used  by  the  stability  code  linear. x  (Herbert  1990)  and 
the  more  advanced  and  efiicient  stability  code  LISA  developed  at  DynaFlow. 
These  files  can  be  directly  incorporated  in  the  PSE  code. 

A  different  approach  has  been  developed  for  the  compressible  case.  The 
full  nonlinear  equations  which  account  for  high  Mach  numbers  (up  to  quartic 
terms)  are  formulated  for  Cartesian  coordinates.  This  challenging  set  of  equa¬ 
tions  is  used  for  all  configurations.  Additional  subroutines  and  insert  files 
have  been  incorporated  to  account  for  the  metric  by  a  first  transformation 
from  the  coordinates  of  the  basic-flow  computation  to  Cartesian  coordinates 
and  a  second  transformation  from  Cartesian  coordinates  to  the  coordinate 
system  for  the  stability  computation.  The  strict  separation  of  the  metrics 
from  the  physics  makes  the  problem  more  manageable.  The  distinction  of  two 
transformations  allows  greater  flexibility  in  adapting  to  the  results  of  CFD 
codes  with  structured  grids  and  choosing  optimal  coordinates  for  the  PSE 
analysis.  The  choice  of  proper  coordinates  for  solving  the  PSE  is  important 
since  the  normal-mode  solutions  of  the  linear  stability  and  the  extended  form 
of  the  PSE  solution  involve  approximations  that  depend  on  the  coordinate 
system. 
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While  concepts,  modular  structure,  and  many  of  the  insert  and  data  files 
are  similar  to  previous  codes,  a  special  version  of  LISA  and  the  compressible 
PSE  code  COPS  have  been  developed  in  cooperation  with  G.  Stuckert  of 
DynaFlow.  The  testing  of  these  codes  has  been  challenging  since  there  are 
rarely  any  reliable  data  to  compare.  To  guarantee  correct  function,  we  have 
devised  a  suite  of  tests  that  solves  the  same  group  of  problems  in  different 
coordinate  systems.  The  results  of  these  tests  are  not  expected  to  completely 
agree  but  are  within  the  known  approximations  involved  in  the  analysis. 

For  flows  in  realistic  geometries,  various  procedures  have  been  coded  and 
compared  to  obtain  the  metric  terms  and  flow  quantities  as  well  as  their 
derivatives  as  accurately  as  possible.  This  task  is  non-trivial  since  the  finite- 
difference  methods  used  for  calculating  the  basic  flow  is  usually  only  of  sec¬ 
ond  order  in  space.  There  are  certain  trade-offs  between  solving  the  stability 
problem  with  spectral  or  finite-difference  methods.  The  spectral  method  is 
advantageous  by  requiring  lower  derivatives  than  the  widely  used  4th-order 
compact  scheme.  The  compact  scheme,  however,  can  utilize  the  given  data  at 
the  grid  points  directly  while  the  spectral  method  uses  less  points  in  a  differ¬ 
ent  yet  rigid  distribution.  Suitable  procedures  have  been  developed  for  both 
cases  since  the  spectral  method  is  far  superior  for  calculating  eigenvalue  spec¬ 
tra  needed  to  identify  the  complete  set  of  unstable  modes.  Finite-difference 
methods,  however,  are  more  efficient  for  tracing  a  specific  mode  over  variable 
parameters  and  for  the  marching  procedure  of  the  PSE  code. 

Formulation  of  the  PSE  equations  and  practical  applications  reveal  two 
severe  disadvantages  of  the  compact  method:  the  need  for  writing  the  equa¬ 
tions  as  a  first-order  system  and  the  additional  differentiation  of  the  equations 
required  to  reduce  the  truncation  error.  The  first-order  system  introduces 
the  inverse  of  the  transverse  velocity  component  of  the  basic  flow  into  the 
coefficients  of  the  differential  equations.  This  transverse  component  is  gen¬ 
erally  small  and  difficult  to  accurately  obtain  by  finite-difference  methods. 
Moreover,  the  transverse  component  can  change  sign  and  cause  undesirable 
singularities  of  the  problem.  We  have  therefore  performed  extensive  studies 
of  high-order  numerical  methods  suitable  for  boundary-layer,  stability,  and 
PSE  calculations. 


5 


2.2  Hermitian  Methods 

The  goal  of  our  efforts  was  to  find  finite  difference  schemes  that  are  of  high 
order,  can  be  applied  to  second-  and  third-order  differential  equations,  lead 
to  block-tridiagonal  systems  (for  efficient  solution),  and  work  for  nonuniform 
grids.  The  result  of  this  work  was  the  selection  of  Hermitian  methods  which 
involve  only  three  neighboring  points.  The  use  of  the  function  values  and 
derivatives  at  these  points  introduces  enough  degrees  of  freedom  to  reduce 
the  truncation  error  to  a  desirable  order. 

Hirsh  (1975)  has  used  Hermitian  high-order  finite-difference  schemes  to 
solve  some  fluid  mechanics  problems.  He  attributes  the  scheme  to  “a  sugges¬ 
tion  made  by  Kreiss,  pertaining  to  a  new  compact  differencing  of  fourth-order 
accuracy.”  The  idea  of  formulating  high-order  finite-difference  approxima¬ 
tions  in  a  compact  way,  i.e.  by  utilizing  the  values  of  derivatives  instead  of  an 
increasing  number  of  neighboring  function  values,  dates  in  fact  back  to  the 
work  of  Hermite  (1912).  Selected  formulas  are  tabulated  by  Collatz  (1966, 
Table  III,  pp.  538ff)  who  also  discusses  Hermite’s  method.  Both  Collatz  and 
Hirsh  give  expressions  only  for  equidistant  grid  points. 

To  maintain  high  order,  special  attention  has  to  be  paid  to  the  boundary 
conditions.  To  achieve  the  same  4th-order  approximation  for  interior  and 
boundary  points,  Hirsh  uses  boundary  formulas  that  involve  more  than  the 
three  neighboring  points  in  the  interior,  although  this  procedure  prevents  the 
use  of  the  highly  efficient  solvers  for  block-tridiagonal  systems. 

The  use  of  Hermitian  methods  on  coarse  and  nonuniform  grids  is  dis¬ 
cussed  by  Adam  (1975).  Adam  uses  boundary  formulas  of  lower  order  than 
in  the  interior  and  uses  nonuniform  grids  to  improve  the  accuracy  near  the 
boundaries.  He  concludes  that  the  lower-order  boundary  formulas  have  no 
detrimental  effect  on  the  accuracy  of  the  solution.  However,  Adam’s  uti¬ 
lization  of  nonuniform  grids  is  in  conflict  with  our  desire  to  use  the  grid 
spacing  to  account,  and  increase  the  resolution,  for  physical  phenomena  such 
as  internal  boundary  layers. 

We  have  developed  a  systematic  procedure  to  generate  all  linearly  inde¬ 
pendent  Hermitian  formulas  using  up  to  the  third  derivative.  For  a  system 
of  N  second-order  differential  equations,  two  of  these  formulas  are  needed  at 
every  interior  point,  three  at  the  boundary.  By  exploiting  the  null  space  of 
the  governing  system  of  Taylor  expansions,  we  have  constructed  methods  of 
at  least  sixth  order  which  increases  to  seventh  order  as  the  number  of  grid 
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points  increases.  By  an  elimination  procedure,  the  original  dimension  3N  of 
the  blocks  in  the  diagonal  system  can  be  reduced  to  2N,  the  same  number 
as  with  the  compact  method.  Unavoidable  deviations  from  tri-diagonality 
at  the  boundaries  require  a  dedicated  subroutine  for  solving  the  algebraic 
system.  In  contrast  to  the  usual  solvers  for  banded  matrices,  our  subroutine 
exploits  the  block  structure  and  increases  the  accuracy  of  the  solution  by 
partial  pivoting  within  blocks. 

The  Hermitian  method  has  been  completely  developed  for  systems  of 
second  and  third  order.  Special  elimination  procedures  have  been  written  for 
systems  derived  from  the  Navier-Stokes  equations  which  have  no  boundary 
conditions  on  the  pressure.  The  method  matches  the  accuracy  of  a  spectral 
method  at  a  fraction,  say  2-5%  of  the  computer  time  and  reduced  demand 
for  memory.  The  method  is  key  to  using  the  PSE  code  on  workstations,  a 
definite  advantage  for  practical  applications.  A  further  advantage  is  the  use 
of  the  governing  equations  in  original  form  without  conversion  to  a  first-order 
system. 

A  report  on  this  method  including  applications  to  challenging  test  cases 
and  convergence  studies  has  been  prepared  and  will  be  refined  for  publication 
as  time  permits. 

2.3  Stability  of  3D  Boundary  Layers 

Detailed  studies  have  been  conducted  to  find  the  correct  solution  to  the 
problem  of  normal-mode  evolution  in  3D  boundary  layers.  Controversial 
views  on  this  subject  have  been  presented  by  Cebeci,  Mack,  Malik,  and 
Nayfeh  in  their  papers  and  at  a  workshop  of  the  US  Transition  Study  Group. 
At  this  time,  a  clear  answer  has  been  found  only  for  boundary  layers  with 
3D  velocity  over  2D  geometry,  e.g.  infinite  or  conical  wings.  Our  findings 
are  in  essential  agreement  with  Mack.  For  an  infinite  wing  the  spanwise 
wavelength  of  disturbances  must  remain  constant.  This  constraint  is  violated 
by  most  stability  codes  for  practical  applications,  in  particular  by  all  optional 
strategies  of  COSAL.  We  plan  to  analyze  the  effect  of  the  2D  approximation 
by  comparative  analysis  of  the  same  boundary  layer  on  an  infinite  and  a 
conical  wing. 

For  fully  3D  boundary  layers,  we  have  performed  various  conceptual  stud¬ 
ies  aiming  at  the  linear  stability  (as  a  first  step)  of  flows  over  highly  tapered 
and  twisted  wings  or  turbine  blades  with  hub  and  tip.  These  studies  indicate 
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that  the  approach  will  be  necessarily  computational  since  the  basic  flow  can 
only  be  provided  by  CFD.  Two  different  computational  methods  have  been 
drafted  which  differ  in  the  spectral  vs.  finite-difference  representation  on 
spanwise  direction.  Implementation  and  validation  of  these  methods  would 
require  at  least  two  man  years  and  is  beyond  the  working  period  and  budget 
for  this  contract. 

2.4  Support  Codes 

An  unexpected  amount  of  work  has  been  required  to  analyze  and  satisfy  the 
prerequisites  for  obtaining  the  basic  flow  for  stability  and  transition  analy¬ 
sis  in  realistic  cases.  The  problem  usually  starts  with  the  geometry  of  the 
body  given  as  an  extensive  set  of  truncated  coordinates.  The  metric  terms 
of  the  body-fitted  coordinate  system  must  be  extracted  from  these  coordi¬ 
nates.  The  dense  spacing  together  with  the  loss  of  significant  digits  causes 
intolerable  errors  in  metric  terms  that  involve  higher  derivatives  of  the  body 
slope  with  respect  to  the  arc  length.  Various  utilities  have  been  written  and 
tested  to  smooth  the  data  without  loosing  significant  information.  While 
these  efforts  were  successful  in  some  cases,  the  procedures  failed  in  other 
situations.  Since  stability  analysis  deals  with  the  sensitivity  of  the  flow  to 
small  disturbances,  significant  results  can  only  be  obtained  by  keeping  the 
input  data,  here  the  geometry,  free  of  unnecessary  disturbances  by  round-off 
errors.  Once  information  is  lost,  it  cannot  be  retrieved  by  densely  spacing  in¬ 
accurate  coordinates.  In  this  case,  a  sparser  set  of  coordinates  appears  more 
appropriate.  In  any  case,  careful  attention  to  details  of  geometry,  curvature, 
and  surface-pressure  distribution  are  key  to  a  successful  stability  analysis. 

The  second  major  problem  area  is  the  use  of  Navier-Stokes  solvers  versus 
a  combination  of  an  Euler  solver  with  a  boundary-layer  code.  The  analy¬ 
sis  of  both  Navier-Stokes  solutions  and  Euler  solutions  for  flows  over  wings 
and  turbine  blades  has  revealed  unexpected  flaws  in  these  solutions.  While 
our  access  to  the  variety  of  codes  is  very  limited,  the  present  status  of  com¬ 
puters  and  computational  methods  appears  to  prevent  sufficiently  accurate 
solutions  for  the  boundary  layer.  In  addition,  the  high  cost  of  current  DNS 
is  unaffordable  in  engineering  practice.  This  view  is  shared  by  various  design 
engineers  from  industrial  R&D  labs. 

The  combination  of  an  Euler  code  with  a  boundary-layer  code  is  currently 
seen  as  the  best  way  to  obtain  the  basic  flow.  Of  the  results  of  an  Euler  code. 
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only  the  surface-pressure  distribution  is  actually  needed.  The  surface  pres¬ 
sure  is  considered  the  easiest-to-get  and  most  accurate  result  of  these  codes 
and  has  been  the  goal  of  much  of  the  CFD  development  in  the  past.  Velocities 
and  stream-line  orientation  at  the  surface  can  be  used  for  cross-checks  with 
the  boundary-layer  calculation.  Good  boundary-layer  codes  are  rare  to  find. 
Most  codes  are  rather  old  and  use  finite-difference  methods  of  only  second 
order.  We  have  analyzed  the  special  version  WING  of  the  Kaups-Cebeci  code 
commercially  used  in  combination  with  COSAL  and  found  unacceptable  in¬ 
accuracies  especially  for  wings  with  small  taper.  The  accuracy  was  improved 
by  recoding  sections  of  the  code.  We  have  designed  a  new  code  that  solves  an 
extended  set  of  flow  problems  with  Hermitian  finite-difference  methods.  Ex¬ 
tensions  concern  primarily  the  boundary  conditions  that  account  for  suction 
and  prescribed  temperature  or  heat  flux  at  the  body.  The  implementation  of 
the  design  is  not  intended  within  this  project  since  we  want  to  dedicate  our 
resources  to  the  fundamental  a.spects  of  transition  prediction  in  practice. 

2.5  Transition  Mechanisms 

With  the  PSE  research  code,  we  have  focused  on  establishing  a  catalogue  of 
transition  mechanisms  originating  from  different  input  models  for  2D  and  3D 
boundary  layers.  The  mechanisms  in  2D  boundary  layers  without  wall  cur¬ 
vature  caused  by  secondary  instability  or  direct  mode  interaction  have  been 
analyzed  with  special  emphasis  on  the  details  of  the  coupling  mechanism. 
Prototype  results  are  discussed  in  the  AGARD  course  notes  on  “Parabolized 
Stability  equations”  in  Appendix  A.  For  flows  with  primary  Gortler  or  cross¬ 
flow  instability  the  transition  picture  is  less  well  understood.  Both  types  of 
instability  lead  to  disturbance  growth  with  larger  rates  than  the  feeble  TS 
instability.  Nonlinear  terms  act  stabilizing  at  large  amplitudes  in  excess  of 
5%,  say.  The  internal  shear  layers  that  are  considered  the  cause  of  inviscid 
secondary  instabilities  form  only  in  highly  nonlinear  stages  of  the  evolution. 
The  secondary  instability  analysis  of  these  nonlinear  stages,  though  straight¬ 
forward  by  accounting  for  a  number  of  harmonics  or  using  two-dimensional 
eigenfunctions,  is  computationally  expensive  and  time-consuming  because 
it  has  to  be  repeated  at  numerous  levels  of  the  nonlinear  evolution.  Some 
results  of  such  analyses  have  been  obtained  by  Malik  at  NASA  Langley. 

Most  of  our  analyses  were  conducted  to  obtain  an  overview  over  the  non¬ 
linear  characteristics  of  Gortler  vortices  and  steady  or  unsteady  cross-flow 
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vortices.  Detailed  calculations  were  performed  for  the  swept-wing  flows  ex¬ 
perimentally  studied  by  Radetztski  et  al.  at  ASU  (Saric  1993).  Although  the 
experimental  flow  is  practically  incompressible,  the  compressible  code  COPS 
was  employed  to  account  for  curvature  effects.  The  detailed  comparison  of 
PSE  results  with  the  experimental  data  is  still  ongoing  (M.  Wang,  forthcom- 
ing  thesis).  The  qualitative  agreement  is  obvious.  The  study  of  interactions 
between  nonlinear  cross-flow  vortices  and  unsteady  disturbances  that- should 
lead  to  secondary  instabilities  and  the  breakdown  of  the  disturbed  flow  re¬ 
quired  some  modifications  of  the  PSE  research  code.  Although  these  changes 
have  been  made,  the  computations  with  the  new  version  were  not  completed 
within  the  working  period  of  this  contract. 

2.6  Disturbance  Environment  and  Receptivity 

Currently,  the  PSE  codes  are  developed  to  solve  the  initial-boundary-value 
problems  of  linear  and  nonlinear  disturbance  evolution  in  quasi-3D  (three 
velocity  components  depending  on  two  spatial  coordinates)  flows  over  simple 
or  realistic  geometries. 

To  apply  these  codes  for  practical  benefit,  the  initial  and  boundary  condi¬ 
tions  must  be  specified  to  characterize  the  proper  disturbance  environment  in 
atmospheric  flight,  different  wind  tunnels,  or  gas  turbines.  Very  limited  mea¬ 
surements  are  available  for  wind  tunnels  (typically  the  turbulence  level)  while 
the  operating  environment  of  airplanes  and  gas  turbines  is  largely  unknown. 
The  work  in  this  area  essentially  involves  the  construction  of  model-data  sets 
as  input  to  the  codes  based  on  mechanisms  and  scales  of  instabilities  and 
comparison  of  the  results  with  observed  transition  points.  Uncertainty  arises 
both  from  insufficient  knowledge  of  the  environmental  disturbances  at  the 
boundary  (vibrations,  roughness,  waviness)  or  in  the  free  stream  and  the  re¬ 
ceptivity  of  the  flow  which  determines  the  associated  disturbance  amplitudes 
inside  the  boundary  layer. 

We  have  studied  different  types  of  model  data.  A  first  type  accounts  only 
for  initial  data  while  the  boundary  conditions  remain  homogeneous.  This 
model  closely  follows  the  approach  in  most  of  the  DNS  work  The  model 
consists  of  a  combination  of  stability  modes  with  proper  frequencies  and 
wavelength  (calculated  using  LISA)  and  empirically  specified  amplitudes. 
Use  of  this  model  in  the  analysis  of  a  group  of  problems  associated  with  the 
flow  in  gas  turbines  has  been  largely  successful  (Herbert  et  al.  1993).  The 
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work  also  revealed  some  shortcomings  of  this  first  type  of  input  model  in  flows 
which  have  an  extended  region  of  accelerated  and  hence  stable  flow,  such  as 
the  suction  side  of  turbine  blades.  In  these  cases,  the  initial  disturbances 
decay  although  the  turbulent  fluctuations  in  the  core  flow  should  maintain 
a  certain  disturbance  level  in  actual  flows.  We  have  studied,  and  in  simple 
flows,  implemented  the  concept  of  area  distributed  receptivity  into  a  second 
type  of  input  model.  PSE  runs  for  different  physical  disturbances  at  the  edge 
or  at  the  body  have  been  performed  to  account  for  different  aspects  of  the 
disturbance  environment.  The  successful  implementation  of  inhomogeneous 
boundary  conditions  in  the  PSE  code  was  enabled  by  the  greater  flexibility 
of  the  Hermitian  methods  and  has  lead  to  receptivity  studies  as  a  broad  area 
of  application  for  the  PSE  codes.  Receptivity,  instabilities,  and  transition 
can  be  covered  in  a  single  continuous  PSE  run. 

In  various  runs,  we  have  replaced  the  generation  of  initial  data  by  inho¬ 
mogeneous  conditions  at  the  initial  position.  This  approach  is  similar  e.g.  to 
the  introduction  of  wave  packets  in  experiments  or  DNS  runs.  A  pulse,  for 
example,  introduces  a  two-dimensional  spectrum  of  disturbance  waves,  all  in 
proper  phase  relation  (a  problem  in  Caster’s  (1975)  analysis),  and  allows  de¬ 
tailed  studies  of  the  transient  evolution  and  ultimate  growth  or  decay.  A  PSE 
run  to  “repeat”  the  experiment  of  Caster  and  Crant  (1975)  is  completely  set 
up.  The  run  has  not  yet  been  executed  since  the  post-processing  of  the  data 
will  be  time-consuming  and  we  currently  have  insufficient  disk  space  to  hold 
the  data  on-line. 

A  more  tractable  problem  is  the  disturbance  field  generated  by  a  har¬ 
monic  point  source.  For  this  case,  experimental  data  of  Cilev  et  al.  (1981), 
Kachanov  (1985),  Kendall  (1993a,b),  and  Watmuff  (1993)  are  available  as 
well  as  linear  stability  calculations  by  Mack  (1984)  and  Balakumar  &  Malik 
(1992).  Linear  PSE  runs  have  been  performed  with  more  than  200  modes  to 
provide  sufficient  resolution  near  the  source.  The  still  ongoing  analysis  of  the 
various  data  is  performed  in  cooperation  with  L.  Mack.  Cood  agreement  has 
been  found  between  theoretical  results  for  the  farfield  provided  the  heuris¬ 
tic  energy  correction  is  suppressed.  This  correction  has  been  introduced  by 
Gaster  to  “account”  for  nonparallelism  and  led  to  a  better  agreement  between 
theoretical  and  experimental  data  in  previous  comparisons.  The  current  pic¬ 
ture  of  the  experimental  data  is  somewhat  confusing.  From  the  few  DNS 
solutions  for  the  near  field,  it  has  become  clear  so  far  that  size  and  type  of 
the  source  have  an  unexpectedly  strong  influence  on  the  field.  These  result 
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are  relevant  for  the  modeling  of  disturbances  such  as  roughness  elements  ex¬ 
posed  to  sound  or  actuators  in  flow  control  applications.  The  DNS  results 
also  show  certain  shortcomings  in  the  initial  data  for  the  PSE  analysis.  We 
expect  that  larger  DNS  runs  for  the  different  experimental  conditions  will 
help  to  clarify  the  picture.  These  runs  will  be  performed  and  analyzed  once 
additional  computer  resources  (memory,  disk  space)  become  available. 

Other  studies  relate  to  the  possible  causes  for  the  characteristic  differences 
between  wave  packets  originating  from  free-stream  turbulence  and  pulsed 
blowing/suction  at  the  wall  (Kendall  1993).  It  is  likely  that  the  narrow  free- 
stream  packets  “ride”  on  a  streaky  structure  that  is  simultaneously  created 
by  a  convected  source.  The  source  can  only  be  explained  as  a  statistical 
event  in  the  free-stream  turbulence.  We  have  performed  some  PSE  analyses 
of  models  for  this  situation  but  were  unable  to  complete  these  studies  within 
the  working  period. 

2.7  Transient  Disturbance  Growth 

It  has  recently  been  found  (Henningson  &  Schmid  1992,  Butler  &  Parrel 
1992,  Trefethen  1992)  that  disturbances  in  shear  flow  can  grow  to  consid¬ 
erable  amplitudes  even  if  they  are  not  associated  with  unstable  eigenval¬ 
ues  of  the  stability  equations.  This  transient  growth  and  ultimate  decay  is 
strongest  for  disturbances  with  spanwise  periodic  streamwise  vorticity.  The 
linear  growth  process  is  characteristic  of  initial-boundary-value  problems  for 
non-selfadjoint  differential  equations.  While  previous  work  considered  the 
temporal  evolution  of  a  given  initial  field  in  parallel  flows,  we  have  studied 
the  spatial  evolution  of  disturbances  introduced  by  inhomogeneous  bound¬ 
ary  conditions  (spanwise  periodic  suction,  waviness,  or  pressure  variation  at 
the  initial  position)  with  the  PSE  technique.  A  first  set  of  results  has  been 
reported  by  Herbert  &:  Lin  1993.  The  analysis  has  provided  a  detailed  ex¬ 
planation  for  the  occurrence  of  steady  and  low-frequency  Klebanoff  modes 
(Klebanoff  1971)  in  the  flat-plate  boundary  layer.  Depending  on  Reynolds 
number,  frequency,  and  spanwise  wave  number,  boundary  disturbances  can 
create  Klebanoff  modes,  TS  waves,  or  fast-traveling  modes  concentrated  at 
the  edge  of  the  boundary  layer,  as  observed  in  Kendall’s  experiments.  By 
accounting  for  streamwise  wall  curvature,  our  results  show  a  stabilizing  effect 
of  convex  curvature  while  sufficiently  strong  concave  curvature  prevents  the 
decay  of  the  modes:  Klebanoff  modes  are  analytically  connected  to  Gortler 
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vortices  and  continue  to  grow  exponentially.  In  Falkner-Skan-Cooke  flow, 
the  introduction  of  sweep  permits  to  show  the  analytical  connection  between 
Klebanoff  modes  and  cross-flow  vortices.  All  these  modes  share  similar  recep¬ 
tivities.  The  pure  effect  of  pressure  gradients  has  not  yet  been  investigated.  A 
vague  connection  to  the  observed  streamwise  vortices  penetrating  the  stagna¬ 
tion  plane  in  the  flow  over  circular  cylinders  (or  blunt  leading  edges)  requires 
further  strengthening. 

We  have  concluded  that  other  phenomena  such  as  the  “instability”  of  pipe 
and  pipe-entrance  flow  and  the  occurrence  of  cross-flow  vortices  in  rotating 
flows  are  closely  associated  with  the  transient  growth  mechanism.  Analysis 
of  these  phenomena  of  scientific  as  well  as  practical  interest  could  not  be 
pursued  with  the  available  manpower. 

Environmental  Disturbances 

Free-Stream  Disturbances 
Boundary  Disturbances 

'  T  , 

I  Receptivity  Mechanisms 


Linear  Mechanisms 


Figure  1:  Pathways  to  Transition  and  Turbulence. 


13 


The  transient  growth  mechanism  has  been  incorporated  into  a  new  over¬ 
view  of  pathways  to  turbulence  shown  in  figure  1.  This  figure  is  the  draft  of 
the  key  figure  for  a  discussion  of  the  transition  issue  by  Morkovin,  Reshotko 
&  Herbert  to  be  submitted  to  Science. 

Various  other  receptivity  studies  are  in  direct  context  with  flow  control 
and  the  analysis  of  actuator  concepts.  The  insight  into  the  transient  flow 
originating  from  an  actuator  is  important  for  the  layout  of  the  experimen¬ 
tal  arrangement  of  the  flow-control  system.  Details  of  these  studies  will  be 
reported  under  contract  F49620-93-1-0135. 

2.8  Adjoint  Linear  Problems 

The  adjoint  problem  to  the  Orr-Sommerfeld  equation  has  often  been  solved 
to  determine  Landau  constants  in  perturbation  series  for  nonlinear  stability 
problems.  The  theory  of  transient  disturbance  growth  in  shear  flows  heav¬ 
ily  rests  on  the  non-selfadjoint  nature  of  the  initial-boundary-value  problem. 
Zhigulev  &  Fedorov  (1987)  used  the  biorthogonal  system  of  eigenfunctions  (of 
the  original  and  adjoint  problem)  to  study  the  boundary-layer  receptivity  to 
acoustic  disturbances.  At  the  APS  meeting  1991,  Hill  reported  on  a  theoret¬ 
ical  approach  for  analyzing  the  restabilization  of  the  wake  behind  a  circular 
cylinder  based  on  the  solution  to  the  adjoint  linearized  Navier-Stokes  equa¬ 
tions  (Hill  1992).  This  latter  work  showed  that  the  adjoint  problem  possesses 
not  only  formal  mathematical  properties  but  also  a  deep  physical  meaning: 
the  adjoint  solution  describes  the  sensitivity  (or  receptivity)  of  the  flow  to 
small  forcing  of  different  nature.  Therefore,  the  adjoint  solution  provides  a 
rigorous  tool  to  the  analysis  of  receptivity  and  flow-control  strategies. 

We  have  started  to  investigate  the  use  of  adjoint  solutions  and  to  formu¬ 
late  and  code  the  equations  in  increasingly  general  situations.  The  initial 
work  is  for  incompressible  3D  boundary  layers.  Starting  from  the  linearized 
Navier-Stokes  equations,  it  becomes  immediately  clear  that  adjoint  solutions 
can  be  generated  for  the  local  analysis  in  parallel  flows  solving  ordinary 
differential  equations.  In  nonparallel  boundary  layers  the  adjoint  PSE  can 
be  formulated  and  solved  by  marching  upstream.  For  given  wave  numbers 
and  frequency  and  the  associated  adjoint  eigenfunction  at  some  initial  posi¬ 
tion,  contour  lines  of  appropriate  components  of  the  solution  directly  indicate 
where  to  place  a  disturbance  generator  to  obtain  the  strongest  disturbance  at 
the  (downstream)  initial  position.  Application  to  the  Blasius  boundary  layer 
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provides,  e.g.  the  intuitively  clear  recommendation  to  place  the  vibrating 
ribbon  near  the  critical  layer  at  branch  I  to  obtain  the  largest  TS  wave.  An¬ 
alyzing  proper  integrals  involving  the  adjoint  solution  provides  the  effects  of 
finite-size  weak  sources  of  momentum,  mass,  or  vorticity  in  the  interior  or  at 
the  boundaries.  The  effects  of  pressure  sources  are  not  yet  fully  understood. 

Hill  (1993)  has  recently  presented  results  for  receptivity  in  parallel  flows 
and  locally  parallel  boundary  layers.  He  has  also  used  the  2D  PSE  to  compute 
the  adjoint  solution  in  the  nonparallel  Blasius  boundary  layer  (Hill  1993b). 
We  consider  a  coordination  of  efforts  or  a  cooperation  on  this  promising  new 
area.  Since  his  continued  funding  at  CTR  is  uncertain  and  receptivity  issues 
are  not  a  focus  of  this  center.  Dr.  Hill  may  be  available  to  work  with  us 
either  at  the  University  or  at  DynaFlow. 

2.9  Flow  Control  and  Neural  Networks 

During  a  part  of  the  working  period,  our  activities  have  been  redirected  from 
the  analysis  of  transition  to  the  control  of  transition  in  boundary  layers. 
The  work  of  two  students  in  this  area  is  supported  by  AFOSR-AASERT 
Fellowships. 

Our  capability  to  simulate  transitional  flows  efficiently  has  been  exploited 
to  prepare  a  joint  theoretical/computational/experimental  study  of  flow  con¬ 
trol.  The  experiments  are  conducted  by  J.  Haritonidis.  Under  this  contract, 
we  have  concentrated  on  neural  networks  for  flow  control,  specifically  on 
suitable  architecture,  training  methods,  fault  tolerance,  and  self-adaptivity. 
Various  nonlinear  dynamical  systems  such  as  the  Lorenz  equations  have  been 
modeled  with  as  simple  as  possible  networks.  The  Neuralworks  Professional 
II  software  has  been  acquired  to  simulate  networks  on  one  of  our  worksta¬ 
tions.  With  three  inputs,  three  outputs,  and  six  “neurons”  in  a  single  hidden 
layer,  the  Lorenz  equations  can  be  well  modeled.  Besides  standard  dynam¬ 
ical  systems,  we  have  successfully  modeled  a  nonlinear  feedback  system  for 
engine  speed  control. 

For  training,  we  had  the  best  experience  with  DRS  (direct  random  search) 
which  is  relatively  slow  but  does  not  settle  in  local  minima  like  back  prop¬ 
agation.  A  problem  that  virtually  found  little  attention  in  the  literature  is 
to  measure  the  quality  of  a  trained  network  or  to  compare  the  quality  of 
different  networks.  The  residual  rms  error  is  not  a  good  quality  indicator. 
We  have  used  various  statistical  methods  to  analyze  the  discrete  response 
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of  the  original  system  and  of  the  neural  network  for  the  same  input  data. 
While  the  development  is  continuing,  we  have  now  tools  at  hand  to  compare 
network  function  and  to  measure  the  deterioration  of  the  network  response 
for  increasing  prediction  times  and  the  degradation  through  failure  of  con¬ 
nections  or  complete  neurons.  Interestingly,  the  loss  of  a  heavily  weighted 
connection  can  be  more  devastating  than  the  loss  of  the  whole  neuron  to 
which  it  is  connected.  A  report  on  these  studies  has  been  given  at  the  AIAA 
Shear-Flow  Conference  in  Orlando  (Fan  &  Herbert  1993). 

After  studies  on  fundamentals,  we  have  concentrated  on  a  flow-control 
situation  that  can  be  both  simulated  and  tested  in  a  wind  tunnel.  For  sim¬ 
plicity,  we  have  chosen  the  2D  problem  of  wave  packets  periodically  generated 
by  a  ribbon  or  point  source.  Each  packet  consists  of  a  mix  of  2D  waves  with 
different  frequencies.  The  mixture  of  phases  and  phase  speeds  poses  a  far 
more  challenging  control  problem  than  the  simple  wave  cancellations  in  ear¬ 
lier  work.  The  test  case  has  been  used  to  study  number  and  arrangements  of 
sensors,  actuators,  and  feedback  connections.  The  preliminary  design  devel¬ 
oped  by  computer  has  been  successfully  tested  in  the  experiment.  In  this  way, 
we  have  not  only  tested  the  concept  and  function  of  neural- network  control 
but  also  our  design  capabilities  prior  to  attacking  more  complex  situations. 

3  Personnel 

The  following  personnel  has  participated  in  the  work  and  has  been  partially 
supported  under  this  contract: 

Th.  Herbert,  Principal  Investigator 

Rihua  Li,  Postdoctoral  Associate 

Mengjie  Wang,  Ph.D.  Student 

Xuetong  Fan,  Ph.D.  Student 

Eugene  Kalinin,  M.S.  Student 

Charlotte  Herbert,  Systems  Programmer  2 

Lorenz  Hofmann  (MS  Student)  has  cooperated  in  this  program  under  an 
AFOSR-AASERT  Fellowship.  Rihua  Li  who  accepted  another  position  in 
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March  1993  and  Mengjie  Wang  have  worked  on  aspects  and  applications  of 
the  PSE.  M.  Wang  will  receive  his  Ph.D.  in  the  near  future.  Xuetong  Fan  and 
Lorenz  Hofmann  have  worked  on  neural  networks  and  their  application  to 
flow  control.  Their  work  continues  under  contract  F49620-93-1-0135.  Eugene 
Kalinin  cooperates  in  the  analysis  and  physical  interpretation  of  solutions  to 
the  adjoint  linearized  Navier-Stokes  equations. 

4  Publications 

The  following  publications  were  completed  or  originated  from  work  under 
support  by  this  contract: 

“Nonlinear  Evolution  of  Secondary  Instabilities  in  Boundary-Layer  Transi¬ 
tion,”  by  J.  D.  Crouch  and  Th.  Herbert,  J.  Theor.  Comp.  Fluid  Dynamics 
(1992),  in  press. 

“Exploring  Transition  by  Computer,”  by  Th.  Herbert,  J.  Appl.  Num. 
Math.,  Vol.  7,  No.  1,  pp.  3-27  (1991). 

“Analysis  of  the  Linear  Stability  of  Compressible  Boundary  Layers  using 
the  PSE,”  by  F.  P.  Bertolotti  and  Th.  Herbert,  J.  Theor.  Comp.  Fluid 
Dynamics,  Vol.  3,  pp.  117-124  (1991). 

“A  Note  on  the  Calculation  of  Landau  Constants,”  by  J.  D.  Crouch  and  Th. 
Herbert,  Phys.  Fluids  A,  Vol.  5,  pp.  283-285,  (1993). 

“Computation  of  Laminar  Flow  over  a  Long  Slender  Axisymmetric  Blunted 
Cone  in  Hypersonic  Flow,”  by  V.  Esfahanian,  Th.  Herbert,  and  0.  R. 
Burggraf,  AIAA  Paper  No.  92-0756  (1992). 

“Stability  of  Hypersonic  Flow  over  a  Blunt  Body,”  by  Th.  Herbert  and 
V.  Esfahanian,  Proc.  AGARD  Symposium  Theoretical  and  Experimental 
Methods  in  Hypersonic  Flows,  Torino,  Italy.  (1992).  To  appear. 

“Linear  and  Nonlinear  Stability  of  the  Blasius  Boundary  Layer,”  by  F.  P. 
Bertolotti,  Th.  Herbert,  and  P.  R.  Spalart,  J.  Fluid  Mech.  Vol.  242,  pp. 
441-474  (1992). 

“Stability  and  Transition  on  Swept  Wings,”  by  G.  K.  Stuckert,  Th.  Herbert, 
and  V.  Esfahanian,  AIAA  Paper  No.  93-0078  (1993). 
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“Effects  of  Free-Steam  Turbulence  on  Boundary-Layer  Transition,”  by  Th. 
Herbert,  G.  K.  Stuckert,  and  V.  Esfahanian,  AIAA  Paper  No.  93-0488 
(1993). 

“Parabolized  Stability  Equations,”  by  Th.  Herbert,  in:  Special  Course  on 
Progress  in  Transition  Modeling,  AGARD  Report  No.  793  (1993). 

“Studies  on  Boundary-Layer  Receptivity  with  Parabolized  Stability  Equa¬ 
tions,”  by  Th.  Herbert  and  Nay  Lin,  AIAA  Paper  No.  93-3053  (1993). 

“Active  Flow  Control  with  Neural  Networks,”  by  X.  Fan,  L.  Hofmann  and 
Th.  Herbert,  AIAA  Paper  No.  93-3273  (1993). 

“Nonlinear  Analysis  of  Swept  Wing  Transitional  Boundary  Layers,”  by  G. 
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1  SUMMARY 

The  parabolized  stability  equations  (PSE)  are  a  new 
approach  to  analyze  the  streamwise  evolution  of  sin¬ 
gle  or  interacting  Fourier  modes  in  weakly  nonparallel 
flows  such  as  boundary  layers.  The  concept  rests  on 
the  decomposition  of  every  mode  into  a  slowly  varying 
amplitude  function  and  a  wave  function  with  slowly 
varying  wave  number.  The  neglect  of  the  small  sec¬ 
ond  derivatives  of  the  slowly  varying  functions  with 
respect  to  the  streamwise  variable  leads  to  an  initial¬ 
boundary- value  problem  that  can  be  solved  by  nu¬ 
merical  marching  procedures.  The  PSE  approach 
is  valid  in  convectively  unstable  flows.  The  equa¬ 
tions  for  a  single  mode  are  closely  related  to  those 
of  the  traditional  eigenvalue  problems  for  linear  sta¬ 
bility  analysis.  However,  the  PSE  approach  does  not 
exploit  the  homogeneity  of  the  problem  and,  there¬ 
fore,  can  be  utilized  to  analyze  forced  modes  and  the 
nonlinear  growth  and  interaction  of  an  initial  distur¬ 
bance  field.  In  contrast  to  the  traditional  patching 
of  local  solutions,  the  PSE  provide  the  spatial  evo¬ 
lution  of  modes  with  proper  account  for  their  his¬ 
tory.  The  PSE  approach  allows  studies  of  secondary 
instabilities  without  the  constraints  of  the  Floquet 
analysis  and  reproduces  the  established  experimental, 
theoretical,  and  computational  benchmark  results  on 
transition  up  to  the  breakdown  stage.  The  method 
matches  or  exceeds  the  demonstrated  capabilities  of 
current  spatial  Navier-Stokes  solvers  at  a  small  frac¬ 
tion  of  their  computational  cost.  Recent  applications 
include  studies  on  localized  or  distributed  receptivity 
and  prediction  of  transition  in  model  environments 
for  realistic  engineering  problems.  The  following  de¬ 
scribes  the  basis,  intricacies,  and  some  applications 
of  the  PSE  methodology. 

2  Introduction 

The  analysis  of  the  transition  process  in  shear  flows 
is  of  both  fundamental  and  practical  interest.  Over 


the  past  decade,  much  insight  into  this  process  has 
been  gained  from  combined  experimental,  theoretical, 
and  computational  studies.  However,  this  progress 
has  little  impact  yet  on  the  engineering  methods 
used  for  aerodynamic  design,  since  the  successful  re¬ 
search  tools  suffer  from  some  shortcomings  that  pre¬ 
vent  their  application  to  practically  relevant  phenom¬ 
ena.  The  traditional  engineering  tools  still  rest  on 
questionable  concepts  and  rely  on  spotty  empirical 
data  bases  which  are  costly  to  generate  for  advanced 
projects. 


Figure  1:  Stability  diagram  for  the  flat-plate  bound¬ 
ary  layer. 

In  low-speed  experiments  in  wind  tunnels  of  very 
low  turbulence  levels,  “natural”,  i.e.  uncontrolled 
and  unavoidable  transition  in  the  flat-plate  bound¬ 
ary  layer  occurs  at  length-Reynolds  numbers  Re^  be¬ 
low  4  •  10®  or  slightly  higher  in  anechoic  tunnels.  A 
glance  at  the  stability  diagram  for  this  flow  in  fig¬ 
ure  1  shows  that  the  occurrence  of  natural  transi¬ 
tion  restricts  experimental  studies  of  the  controlled 
disturbance  evolution  from  branch  I  to  branch  II  to 
Reynolds  numbers  Re  =  Rex  well  below  Re  <  2000 
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and  hence  to  frequencies  F  >  25.  In  fact,  the  micro¬ 
scopic  experiments  are  restricted  to  this  range,  where 
two-dimensional  waves  are  “most  dangerous”  [1],  On 
the  other  hand,  the  common  criterion  for  the  am¬ 
plitude  ratio  and  the  correlation  with  experimental 
data  show  transition  typically  at  >  10  [2,  3].  To 
obtain  an  amplitude  ratio  of  at  branch  II,  the  fre¬ 
quency  must  be  F  <  25.  At  these  lower  frequencies, 
two-dimensional  waves  are  no  longer  preferred.  Con¬ 
sequently,  the  microscopic  experiments  do  not  nec¬ 
essarily  reveal  the  transition  process  as  it  occurs  in 
practice,  let  alone  the  different  disturbance  environ¬ 
ments  in  free  flight  and  in  wind  tunnels. 

Theoretical  and  computational  studies  of  transi¬ 
tion  and  its  inherent  mechanisms  have  led  to  con¬ 
siderable  progress  in  explaining  and  reproducing  the 
results  of  the  ribbon-controlled  experiments.  Accord¬ 
ingly,  these  studies  have  revolved  around  the  same  pa¬ 
rameter  range  of  relatively  small  Reynolds  numbers, 
high  frequencies,  and  low  growth  rates.  The  linear 
theory  of  primary  instability  is  not  restricted  to  this 
range,  and  the  underlying  parallel-flow  approxima¬ 
tion  is  commonly  assumed  to  be  better  justified  at 
lower  frequencies  and  the  associated  higher  Reynolds 
numbers  downstream  of  branch  L  Therefore,  the 
calculations  are  performed  without  major  concern. 
However,  nonparallel  effects  are  indeed  significant 
for  certain  three-dimensional  disturbances,  invalidate 
Squire’s  theorem  and  the  normal-mode  concept  for 
long-wave  disturbances,  and  may  accumulate  in  the 
N  factor  integration.  Even  if  stream  wise  variations 
of  the  boundary  layer  are  properly  included,  the  pri¬ 
mary  instability  is  only  one  of  the  important  steps 
toward  transition. 

The  development  of  the  theory  of  secondary  insta¬ 
bility  [4]  has  advanced  our  insight  into  the  next  step 
where  strong  vortical  mechanisms  are  activated.  In 
spite  of  some  necessary  approximations,  the  theory 
provides  results  in  good  agreement  with  experiments 
and  enables  systematic  studies  in  a  multi-dimensional 
parameter  space.  However,  the  potential  for  signifi¬ 
cant  improvements  in  predicting  the  onset  of  transi¬ 
tion  is  disappointing  since  the  theory  is  still  linear  and 
hence  cannot  reveal  the  nonlinear  coupling  of  primary 
and  secondary  disturbances  that  appears  crucial  for 
transition.  A  perturbation  scheme  has  been  devel¬ 
oped  to  successfully  analyze  this  nonlinear  stage  [5], 
but  the  method  is  not  as  well  suited  for  routine  appli¬ 
cations  as  many  other  weakly  nonlinear  models  and 
methods  in  this  held. 

Over  the  past  decade,  details  and  structure  of  the 
transition  process  have  been  studied  with  asymptotic 
theories  for  high  Reynolds  numbers  [6].  These  theo¬ 
ries  enlist  impressive  results  for  some  linear  and  non¬ 
linear  problems  and  local  phenomena  that  can  be  cap¬ 
tured  within  a  single  asymptotic  structure.  The  dras¬ 
tically  simplified  equations  often  allow  closed-form  so¬ 
lutions  or  ease  numerical  solutions  to  nonlinear  prob¬ 


lems.  However,  the  quantitative  features  of  the  tran¬ 
sition  process  that  critically  depend  on  the  physical 
parameters,  especially  the  Reynolds  number,  may  not 
be  accessible  to  these  asymptotic  methods.  As  in 
other  asymptotic  theories,  it  cannot  be  a  priori  as¬ 
sured  that  expansion  parameters  such  as  e  — 
are  sufficiently  small  for  quantitative  conclusions. 

With  increasing  capability  and  availability  of  su¬ 
percomputers,  efforts  have  broadened  to  analyze  the 
transition  process  by  direct  numerical  solution  of  the 
Navier-Stokes  equations  (DNS)  [7].  Computer  sim¬ 
ulations  provide  a  wealth  of  detailed  data,  resolved 
in  3D  space  and  time,  for  the  idealized  prototype 
problems  of  traditional  stability  theory,  and  can  be 
viewed  as  a  new  type  of  “experiment”  under  well- 
controlled  conditions.  Extracting  insight  requires  te¬ 
dious  post-processing  yet  is  facilitated  by  the  easy 
access  to  accurate  information  on  quantities  that  are 
difficult  to  measure,  e.g.  vorticity.  The  effects  of  ini¬ 
tial  conditions  or  changes  in  geometry  and  basic  flow 
on  the  transition  location  would  be  easy  to  extract. 
However,  DNS  is  not  yet  feasible  as  an  engineering 
method.  The  enormous  computational  costs  in  both 
memory  and  time  prohibit  parameter  studies  and 
have  constrained  transition  simulations  in  boundary 
layers  to  a  few  runs,  typically  for  the  conditions  of  the 
benchmark  experiments  [8,  9,  10].  All  simulations  are 
for  simple  geometries,  and  only  a  few  track  the  spa¬ 
tial  evolution  in  streamwise  varying  boundary  layers. 
Although  current  research  focuses  more  on  realistic 
simulations  and  algorithmic  improvements,  both  re¬ 
search  and  engineering  practices  demand  less  expen¬ 
sive  methods  without  intolerably  compromising  the 
quality  of  results.  This  demand  is  particularly  strong 
for  compressible  flows,  where  the  computational  re¬ 
quirements  increase  easily  by  an  order  of  magnitude. 
Even  the  Floquet  analysis  of  secondary  instability  in 
a  compressible  flow  can  be  a  challenge  for  a  super¬ 
computer. 

Here  we  discuss  a  less  expensive  method  that  has 
shown  great  potential  for  efficient  transition  simu¬ 
lations.  This  method  rests  on  parabolized  stability 
equations  (PSE)  and  takes  advantage  of  the  efficient 
methods  available  for  solving  parabolic  partial  differ¬ 
ential  equations  with  marching  procedures.  The  PSE 
approach  can  be  viewed  on  one  hand  as  a  WKBJ 
extension  of  the  theoretical  approach  to  nonparal¬ 
lel  flows  and  on  the  other  cLs  a  simplification  of  the 
Navier-Stokes  equations.  The  approach  is  based  on 
earlier  ideas  for  analyzing  the  stability  of  weakly  non¬ 
parallel  flows  [11]  and  on  the  concept  of  convective 
instability  [12]  which  excludes  significant  upstream 
influence  of  disturbances.  While  originally  intended 
to  capture  nonparallel  effects  on  linear  stability  char¬ 
acteristics  of  single  modes,  the  nonlinear  version  of 
the  PSE  turned  out  to  provide  results  on  secondary 
instability  and  mode  interactions  that  suggest  the  use 
of  this  approach  for  analyzing  transition  up  to  the 
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spike  stage  and  beyond.  The  efficiency  of  the  ap¬ 
proach  and  the  good  agreement  of  PSE  and  DNS 
solutions  encourage  further  development  both  as  a 
research  method  and  as  an  engineering  tool  for  tran¬ 
sition  analysis  and  prediction. 

The  following  describes  the  underlying  formalisms, 
mathematical  aspects,  numerical  procedures,  and 
some  applications  to  problems  of  increasing  complex¬ 
ity  in  incompressible  and  compressible  2D  and  3D 
boundary-layer  flows.  The  applications  to  similar¬ 
ity  solutions  are  chosen  to  exhibit  the  essential  steps 
of  the  analysis  as  well  as  the  potential  of  the  ap¬ 
proach.  Some  applications  to  realistic  geometries  are 
discussed  to  demonstrate  the  feasibility  of  PSE-based 
engineering  methods  and  to  shed  some  light  on  the 
requirements  of  transition  predictions  in  practice. 


3  PSE  for  Incompressible  Flow 


The  motion  of  an  incompressible  fluid  with  constant 
density  p  and  viscosity  p  is  governed  by  the  continuity 
equation 

Vv  =  0  (1) 

and  the  nonlinear  momentum  equations 


^  +  (vV)v  =  -ivp  +  i/V2v,  (2) 

01  p 


where  v  is  the  velocity  vector,  p  the  pressure,  and  p  = 
p/p  the  kinematic  viscosity.  In  Cartesian  coordinates 
X,  y,  z,  we  use  the  velocity  components  v  =  (w,  v,  w) 
s^nd  ^  ^  ^ 


where  i,j,k  are  the  unit  vectors  in  the  direction  of 
respectively. 

We  consider  the  stability  of  steady  laminar 
boundary-layer  flows  over  the  plane  y  =  0  with  edge 
velocities  Ueix)  and  We{x)  in  the  x  and  z  direction, 
respectively.  As  common  for  stability  analysis,  we 
decompose  the  total  flow  field  v,p  into  the  steady 
laminar  basic  flow  V,  P  and  the  deviation  (or  distur¬ 
bance)  v',p': 


V  =  V  -h  v'  ,  p  =  P  -f  p'  .  (4) 

Owing  to  the  uniformity  in  z,  the  steady  basic 
flow  takes  the  form  V  =  {U{Xjy),V{x,y),W{x,y)) 
where  V,  Ux,  and  are  small  and  their  deriva¬ 
tives  UxxiVxy^xx  with  respect  to  x  are  of  order 
0{Re~^)  and  negligible.  Here,  the  Reynolds  number 
is  Rex  =  Urx/u  with  a  proper  reference  velocity  Ur> 
Under  the  boundary-layer  approximation,  the  basic 
flow  is  governed  by 


o 

11 

I 

1 

(5c) 

UWx  +  VWy  =  UWyy  , 

(5d) 

where  Px  =  p{U edU e / dx -{■  WedW^/dx)  is  given  from 
potential-flow  theory.  Equations  (5)  are  parabolic 
partial  differential  equations  and  can  be  solved  as  an 
initial-boundary- value  problem  in  x.  Boundary  con¬ 
ditions  are  the  no-slip  conditions 


ff=K  =  W  =  Oaty  =  0  (6a) 

and  the  asymptotic  matching  with  the  edge  condi¬ 
tions, 

W^Weasy-^oo.  (66) 

The  transverse  velocity  V  is  small  but  in  general 
nonzero  at  the  edge  of  the  boundary  layer.  Of  partic¬ 
ular  interest  for  theoretical  studies  are  the  exact  sim¬ 
ilarity  solutions  to  eqs.  (5)  and  (6).  These  solutions 
are  obtained  from  ordinary  differential  equations  and 
thus  do  not  require  the  specification  of  initial  con¬ 
ditions.  Similarity  solutions  are  the  two-dimensional 
Blasius,  Falkner-Skan,  or  Hiemenz  flows  or  the  three- 
dimensional  Falk ner-Skan- Cooke  or  swept  Hiemenz 
flows.  For  three-dimensional  flows,  the  constant  vis¬ 
cosity  permits  solving  first  the  two-dimensional  prob¬ 
lem  (5a)  to  (5c)  for  U  and  V  and  afterwards  eq.  (5d) 
for  W .  Little  modification  is  required  to  extend  the 
basic  flows  to  two-dimensional  mixing  layers,  jets,  or 
wakes,  or  to  axisymmetric  cases.  Nonsimilar  solu¬ 
tions  can  be  generated  by  numerical  techniques. 

For  the  disturbances  v',  p'  we  obtain  the  “stability 
equations” 

<  +  t;'+u)'=0,  (7a) 

u\  -\-Uu'^+  u'Ux  -\-Vu'y+  v'Uy  +  Wu'^  +  ^p'j. 

+  +  =  -{u'u'^  +  v'u'y  +  w'u'J  ,  (76) 

v[  +  Uv'^  +  Vv'y  +  v'Vy+  Wv',  +  ip; 

+  Vyy  +  V'^^)  =  +  V'v'y  +  W'^^)  ,  (7c) 

w[  +  Uw'^  +  u'Wx  +  Vw’y  -^v'Wy^-  Ww',  +  ip; 

-U{w'^^  +  Wyy+W^^)  =  -{u'  W'^  +  V'  W'y-\-W'  W',)  ,  {1(1) 

where  we  have  neglected  the  small  term  u'\4  in 
(7c).  For  a  given  basic  flow,  and  within  the  boundary- 
layer  approximation,  this  system  is  equivalent  to  the 
N avier-Stokes  equations  and  hence  of  elliptic  type. 

For  the  study  of  instabilities,  we  usually  apply  ho¬ 
mogeneous  boundary  conditions 

u'  =  0,  u'  =  0,  la'  =  0  at  y  =  0  ,  (8a) 


+  V',  =  0  ,  (5a) 


0,  v'  0,  0  as  y  — >  oo 


m 


UUx  +  VUy  =  --Px  +  UUyy  ,  (56) 


since  any  inhomogeneous  conditions  are  already  sat¬ 
isfied  by  the  basic  flow.  In  nonlinear  problems, 
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however,  the  fluctuating  components  of  the  distur¬ 
bance  also  create  a  distortion  of  the  steady  mean  flow 
(u,i5,u;).  Since  the  mean  flow  and  its  distortion  are 
governed  by  the  boundary-layer  equations,  we  must 
admit  v  ^  0  as  y  ^  oo  with 

ziy  0  as  y  —  00  (8c) 

and  =  0.  In  general,  the  nonlinear  stability  equa¬ 
tions  (7)  can  also  be  solved  under  inhomogeneous 
boundary  conditions  to  account  for  disturbances  at 
the  wall  or  in  the  free  stream.  Numerical  simula¬ 
tions  often  use  inhomogeneous  boundary  conditions 
to  simulate  disturbance  generators  such  as  a  vibrat¬ 
ing  ribbon.  Inhomogeneous  boundary  conditions  are 
key  to  the  study  of  forced  oscillations  and  receptivity 
problems. 

In  the  following  we  consider  the  equations  in  nondi- 
mensional  form.  For  convenience  we  maintain  the 
notation  yet  delete  the  factors  1/p  and  replace  u  by 
l/Re,  where  Re  is  formed  with  a  length  scale  6^  re¬ 
lated  to  the  boundary-layer  thickness  and  a  charac¬ 
teristic  velocity  Uq.  Henceforth  the  asterisk  denotes 
dimensional  quantities.  Specific  choices  will  depend 
on  the  application. 

3.1  Linear  Stability  of  Parallel  Flow 

To  illustrate  the  conceptual  difference  between  the 
traditional  local  stability  analysis  and  the  PSE  ap¬ 
proach,  we  consider  the  linear  stability  of  a  parallel 
flow  V  =  (£/(y),0,P7(y)). 

3.1.1  Local  Analysis 

For  a  strictly  parallel  flow,  the  linearized  equations 
(7)  reduce  to  the  linear  stability  equations 


=  0  , 

(9a) 

+  v'Uy  +  wu',  +  p; 

“  **yy  —  0  1 

(9*) 

v't  +  Uv'^  +  Wv',+p'^ 

(9c) 

w[  +  Uw'^  +v'[Vy  +  Ww',  +  p'. 

(9d) 

together  with  the  homogeneous  boundary  conditions 
(8a),  (8b). 

Since  the  coefficients  of  these  differential  equations 
are  either  constant  or  depend  only  on  the  basic  flow 
V(y),  we  can  assume  solutions  in  the  form  of  normal 
modes: 

q'(r,y,z,t)  =  q(y)e‘“"^‘^^-‘^ 


where  q'  =  [u' ,  v' ,  w\  p']'^  is  the  vector  of  flow  vari¬ 
ables  and  q  =  [u,  i-,  w,  p]^  is  the  vector  of  the  asso¬ 
ciated  amplitude  functions.  The  exponential  factor 
with  wavenumbers  a  and  0  in  x  and  z,  respectively, 
and  frequency  u;  describes  the  wave  nature  of  the  so¬ 
lution.  Introducing  the  normal  modes  (10)  into  equa¬ 
tions  (9)  yields  the  normal-mode  equations 


iau  4-  i0w  4-  Dv  =  0  , 

(11a) 

i{aU  4-  0W  —  (jj)u  -f  {DU)v  4-  tap 

(116) 

i{aU  4-  0W  —  Lij)v  4-  Dp 

-^(D^ =0, 

(11c) 

i{QU  +  0W  -uj)w  +  {DW)v  +  i7?p 

(lid) 

with  the  boundary  conditions 

u  =  0,  z;  =  0,  iL>=:0aty  =  0, 

(12a) 

u— *>0,  u— *•0,  u;— ►Oasy  —  oo, 

(126) 

and  D  =  d/dy.  These  equations  represent  a  homo¬ 
geneous  sixth-order  system  of  ordinary  differential 
equations  for  the  natural  variables  u,  v,  w,  and  p. 
Depending  on  the  specific  problem  and  the  numer¬ 
ical  method  of  solution,  various  equivalent  forms  of 
the  normal-mode  equations  are  in  use  [13].  Here  we 
mention  only  the  system 

|( - -  i(aU  +  0W  -  w)](D'  -  Q 

'  -  0-) 

+i{aD^U  +  0D‘^W)}  v  =  0  , 

(13a) 

D^  ^  a^  —  0^ 

Re  ^{otU 0W  —  w)]{aw  — 

0u) 

=  -{aDW  -  0DU)v  , 

(136) 

i(au  4-  0w)  =  —Dv  . 

(13c) 

The  fourth-order  differential  equation  ( 13a)  for  v  with 
homogeneous  boundary  conditions  on  v  and  Dv  is 
the  generalized  form  of  the  Orr-Sommerfeld  problem 
originally  derived  with  W  =  0.  This  homogeneous 
problem  provides  eigenvalues  and  eigensolutions  for 
V.  The  second-order  differential  equation  (13b)  for 
aw  -  0u  with  homogeneous  boundary  conditions  on 
aw  —  0u  is  the  generalized  Squire  problem  that  is  in¬ 
homogeneous  for  given  v  ^  0  but  supports  a  separate 
set  of  eigensolutions  for  v  =  0.  Equation  (13c)  per¬ 
mits  the  calculation  of  u  and  w  once  v  and  aw  —  0u 
are  known. 

The  derivation  of  eqs.  (13)  from  eqs.  (9)  first  in¬ 
volves  the  elimination  of  the  pressure  by  taking  the 


(10) 
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curl  Vx  of  the  momentum  equations  which  pro¬ 
vides  the  vorticity  transport  equations.  The  Orr- 
Sommerfeld  equation  is  obtained  by  subtracting  the 
2:  derivative  of  the  x-vorticity  equation  from  the  x 
derivative  of  the  z-vorticity  equation  and  eliminat¬ 
ing  Wz  via  continuity.  The  Squire  equation  is  the 
z  derivative  of  the  y-vorticity  equation.  Equations 
(13)  offer  various  benefits  for  theoretical  studies,  the 
first  of  which  is  the  separation  of  the  eigenvalue  spec¬ 
trum  into  Orr-Sommerfeld  modes  and  Squire  modes 
that  are  associated  with  different  physical  mecha¬ 
nisms.  Squire  modes  are  stable  and  usually  disre¬ 
garded.  Since  the  interest  in  solving  stability  prob¬ 
lems  usually  is  in  growth  rates,  not  in  the  structure 
of  the  eigenmodes,  the  Squire  equation  (13b)  is  rarely 
used.  In  addition,  the  Orr-Sommerfeld  problem  for 
oblique  waves  with  /?  ^  0  can  be  reduced  to  the  prob¬ 
lem  for  two-dimensional  waves  with  /?  =  0  by  using 
Squire’s  transformation  [13].  Otherwise,  the  prefer¬ 
ence  for  a  specific  formulation  is  largely  determined 
by  the  choice  of  a  numerical  method  for  solution. 
Within  numerical  errors,  the  results  are  independent 
of  the  formulation. 

It  is  worth  noting  that  the  Orr-Sommerfeld  equa¬ 
tion  requires  D^U  and  D^W  which  are  absent  from 
eqs.  (9).  While  this  requirement  is  easy  to  satisfy  in 
studies  of  similarity  solutions,  higher  derivatives  of 
flow  quantities  may  be  inaccurate  or  difficult  to  ob¬ 
tain  if  the  basic  flow  is  computationally  generated. 
Also,  the  separation  of  modes  is  specific  to  incom¬ 
pressible  problems.  In  compressible  problems,  the 
choice  is  between  the  formulation  in  natural  variables 
or  a  first-order  system  of  equations.  In  both  cases, 
the  separation  of  classes  of  modes  requires  a  posieri- 
ori  study. 

Independent  of  the  formulation  and  for  given  U  and 
W,  the  eigenvalue  problem  can  be  reduced  to  solving 
a  complex  characteristic  equation  of  the  form 

=  Q  (14) 

which  yields  two  quantities  per  eigenmode,  provided 
all  others  are  given.  For  boundary  layers,  it  is  most 
appropriate  to  specify  the  real  quantities  i?e,  and 
u  and  determine  a  =  otr^ioci  from  eq.  (14)  to  obtain 
the  spatial  growth  rate  —oti  in  x-direction.  Unfortu¬ 
nately,  different  powers  of  a.  appear  in  the  equations 
which  causes  difficulty  for  many  eigenvalue  solvers. 
Therefore,  the  problem  is  often  solved  with  the  eigen¬ 
value  u  for  real  /2e,  a,  and  P  to  obtain  the  temporal 
growth  rate  w,*. 

Let  us  assume  we  solve  the  homogeneous 
boundary-value  problem  with  Newton’s  method  for 
a  single  eigenmode.  Besides  i?e,  /?,  a;,  /7,  and  W , 
we  provide  an  initial  guess  a  for  the  eigenvalue  and 
q  for  the  eigenfunction,  e.g.  the  results  from  a  previ¬ 
ous  solution  for  different  parameters.  Typically,  the 
Newton  iteration  converges  after  a  few  steps  within 
a  given  convergence  criterion  to  an  approximation 


=  a(i^e, /?,  cj),  q(y)  of  the  eigensolution.  In  a 
strictly  parallel  flow,  this  solution  is  valid  for  any  x. 

In  a  boundary  layer,  the  situation  is  more  compli¬ 
cated.  Owing  to  the  stream  wise  variation  of  the  flow, 
the  stability  characteristics  change,  and  the  compu¬ 
tation  has  to  be  repeated  at  numerous  streamwise 
positions  Xn  for  appropriate  input  data  and  local  ve¬ 
locity  profiles  V.  The  choice  of  wavenumbers  and  fre¬ 
quency  at  a  new  position  is  ambiguous.  It  is  tempting 
to  use  the  parameters  of  the  previous  position,  but 
since  the  length  and  time  scale  usually  involve  the 
x-dependent  boundary-layer  thickness,  a  fixed  value 
of  /?,  for  example,  describes  different  physical  wave¬ 
lengths  at  different  positions.  On  the  other  hand,  the 
stability  analysis  aims  at  finding  the  total  amplifica¬ 
tion  and  the  amplitude  growth  of  the  “most  danger¬ 
ous”  modes  that  evolve  in  the  streamwise  direction 
and  thus  requires  maintaining  the  physical  proper¬ 
ties  of  the  modes.  In  the  traditional  linear  stability 
analysis,  the  physical  connection  between  the  solu¬ 
tion  at  different  streamwise  stations  is  obviously  lost 
by  the  assumption  of  a  parallel  flow  which  implies  the 
same  solution  for  all  x.  Retrieving  this  physical  con¬ 
nection  in  general  three-dimensional  boundary  layers 
is  a  rather  intricate  and  controversial  matter. 

For  the  z-independent  flows  considered  here,  it  can 
be  shown  [14]  that  the  dimensional  frequency  /’*'  and 
spanwise  wavelength  A*  of  a  mode  must  be  indepen¬ 
dent  of  X.  With  these  physical  constraints,  we  can 
calculate  a  sequence  of  eigenvalues  a(x„)  and  obtain 
the  amplitude  ratio  A{x)/Aq  from 

=  (15) 

where  the  asterisk  denotes  dimensional  quantities  and 
Aj  IS  the  amplitude  of  the  mode  at  the  first  neutral 
point  x*f.  Within  a  linear  framework,  the  modes  are 
independent  and  the  “most  dangerous”  modes,  i.e. 
those  with  the  largest  amplitude  ratios  at  any  posi¬ 
tion  can  be  identified  from  the  envelope  of  all  am¬ 
plitude  growth  curves  for  different  values  of  /*  and 

3.2  PSE  Concept  in  Natural  Variables 

The  PSE  concept  was  originally  pursued  [15]  to 
clarify  the  different  results  of  Bouthier  [16,  17], 
Gaster  [11],  and  Saric  &  Nayfeh  [18]  on  the  nonparal¬ 
lel  stability  of  the  Blasius  flow  to  two-dimensional  TS 
waves.  For  nonparallel  flows,  the  normal-mode  solu¬ 
tions  (10)  are  not  strictly  valid  since  the  coefficients 
of  the  linearized  equations  (7)  depend  weakly  on  x. 
Both  the  amplitude  function  q  and  the  wavenumber 
a  change  according  to  nonparallel  effects.  Following 
a  WKBJ  analysis,  we  maintain  the  decomposition  of 
q'  into  an  amplitude  function  q  and  a  wave  function 
X  and  write  the  disturbances  in  the  form 

q'(x,  y,  z,  i)  =  q(x,  2/)x(x,  <)  , 


(16a) 
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X  -  exp(26>(x)  +  i(5z  -  iujt)  ,  (9^  ^  a{x)  .  (166) 

Similar  expressions  have  been  applied  in  the  cited  ear¬ 
lier  work.  We  note  that 

q'r  =  -f-  q^)x  ,  (17a) 

q^rr  -  (-^^q  +  2iaq:^  +  -h  2a3;q)x  ,  (176) 

and  by  substituting  the  components  of  (16)  into  the 
linearized  eqs.  (7)  we  obtain  in  nondimensional  form 

iau  -i-  Dv  +  iPw  +  =  0  ,  (18a) 

\-^[D-  -  [a-  +  /?“)]«  +  i{aU  +0W  -  w)u 


It  is  convenient  to  rewrite  the  parabolized  eqs.  (18) 
in  the  more  compact  form 


Lq  +  +  ^Nq  =  0  , 

ax  ax 


where  the  operators  L,  M,  and  N  act  only  in 
y.  Within  the  PSE  approximation,  the  streamwise 
derivatives  of  q'  take  the  form 


iap  >  -h  {UxU  +  V Du} 


M  =  ,  N 

OCX 


-  (a^  +  j3'^)]v  +  i{aU  +  jSW  -  w)v 


+Dp  ^  +  {DVv  +  VDv] 


2ia 


If  ,  . 


-  {a-  +  p'^)]w  +  i{aU  +  (3W  -  w)r 


These  relations  simplify  the  derivation  of  the  PSE  or 
can  serve  as  a  useful  check. 

Obviously,  we  can  rewrite  eqs.  (18)  in  the  form 
of  a  first-order  system  or  in  a  form  similar  to  the 
Orr-Sommerfeld-Squire  formulation  (12)  [19].  After 
parabolization,  we  obtain  equations  in  the  generic 
form  (19),  (21)  with  different  operators  L,  M,  N . 
While  the  parabolized  first-order  system  involves  ex¬ 
actly  the  same  approximation  as  the  parabolized  eqs. 
(18),  the  additional  differentiations  with  respect  to  x 
necessary  to  obtain  the  Orr-Sommerfeld-Squire  form 
lead  to  a  slightly  different  effect  of  the  PSE  approxi¬ 
mation. 


+DWv  +  iPp}  +  {M4,u  +  VDw) 


+ 


,  (18(f) 


where  now  D  -  d/dy.  The  left-hand  side  of  eqs. 
(18b)  to  (18d)  each  contains  four  groups  of  terms  en¬ 
closed  in  braces.  The  terms  in  the  first  group  are  re¬ 
tained  by  the  normal-mode  equations  (11)  for  parallel 
flow.  The  second  group  originates  from  the  stream- 
wise  changes  of  the  basic  flow:  the  streamwise  gra¬ 
dients  of  U  and  W  and  the  transverse  velocity  com¬ 
ponent  V  associated  with  the  growth  of  the  bound¬ 
ary  layer.  The  third  group  contains  the  streamwise 
derivatives  of  the  amplitude  functions  {u,v^w,p). 
One  of  these  derivatives  also  appears  in  the  conti¬ 
nuity  equation  (18a).  The  last  group  originates  from 
the  streamwise  changes  of  the  wavenumber  a. 

For  a  proper  function  a(x),  the  streamwise  varia¬ 
tion  of  q'  is  governed  by  the  wave  function  x  while  the 
derivatives  q^;  and  are  small.  The  PSE  approxima¬ 
tion  assumes  the  variation  of  q  and  a  as  sufflciently 
small  to  neglect  q^r^,  o^xxi  products  a^q^,  and  their 
higher  derivatives  with  respect  to  x*.  Hence  we  obtain 
the  parabolized  stability  equations  by  neglecting  the 
terms  on  the  right  hand  side  of  eqs.  (18b)  to  (18d). 


3.2.1  Mathematical  Character  of  the  PSE 

At  first  sight,  eqs.  (19)  appear  suitable  for  solution 
as  an  initial-boundary-value  problem.  In  fact,  the 
parabolic  character  of  closely  related  equations  has 
been  mentioned  by  Caster  [11].  Disturbing  is  the 
appearance  of  since  a  affects  the  operators  yet 
is  a  priori  unknown.  The  next  section  will  discuss 
how  to  overcome  this  problem.  The  ad  hoc  numerical 
solution  of  the  PSE  (19)  for  various  two-dimensional 
boundary  layers  is  in  fact  successful  [19].  However, 
numerical  instabilities  occur  for  small  marching  steps 
and  indicate  a  residual  ellipticity  of  the  PSE. 

Ellipticity  of  the  stability  equations  is  necessary 
to  describe  instabilities  in  the  form  of  waves,  e.g. 
TS  waves.  Any  attempt  to  obtain  TS  waves  on 
the  basis  of  the  boundary-layer  equations  or  parab¬ 
olized  Navier-Stokes  equations  (PNS)  must  fail  since 
q^rr  —  is  not  negligible.  The  traditional 

normal-mode  analysis  for  parallel  flows  fully  accounts 
for  the  ellipticity  although  the  separation  of  variables 
leads  to  ordinary  differential  equations.  Through  the 
decomposition  (16),  the  PSE  account  for  the  essential 
part  of  ellipticity  and  remove  it  from  the  equations 
for  the  amplitude  functions  q.  Any  residual  ellipticity 
must  be  associated  with  the  small  deviation  of  waves 
in  nonparallel  flows  from  normal  modes. 
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Depending  on  the  formulation,  different  procedures 
reveal  the  mathematical  character  of  the  partial  dif¬ 
ferential  equations  (19).  We  use  a  modification  of 
the  analysis  presented  by  Haj- Hariri  [20]  for  a  strictly 
parallel  flow  that  is  directly  applicable  to  eqs.  (18). 
We  collect  all  derivatives  of  q  on  the  left  hand  side 
and  indicate  all  other  terms  by  the  dots  on  the  right 
hand  side: 

ciUj;  +  —  0  ,  (22a) 

-^n„  +  Vu,  +  (U--^)u. 

+C2Px  -  =  •••  ,  (226) 

1  2 

H”Py  ~  ,  (22c) 

itc 

1  rr  /rr  ^Za  . 

_—^^^  +  VWy+iU-—)Wx 

5  {22d) 


where  the  constants  Ci  =  l,z  =  0,1,2,  have  been 
introduced  for  convenience.  We  replace  by  Aq  and 
qy  by  /zq  to  obtain  the  system  in  the  form  Aq  =  ... 
with 


(ciA  fi  0  0  ^ 

C  0  0  C2A 

0  C  0  /X 

0  0  C  0  / 


(23a) 


and 


±\/{ReU  —  2za)2  —  —  ReV)]  (26) 

which  are  in  general  complex.  For  the  parabolized 
equations,  co  =  0,  these  roots  reduce  to 


^  _  tijfi  -  ReV) 
ReU  —  2ia 


(27) 


which  for  ar  ^  0  are  still  complex,  though  with  a 
very  small  imaginary  part.  It  is  interesting  to  note 
that  the  sign  of  U  does  not  affect  the  imaginary  parts 
of  the  roots  (25),  (27)  of  the  parabolized  equations. 

The  situation  we  find  for  the  PSE  is  in  fact  very 
similar  to  the  parabolized  Navier-Stokes  equations 
(PNS).  The  equations  can  be  solved  with  march¬ 
ing  techniques  as  an  ill-posed  initial-boundary-value 
problem  provided  the  step  size  is  sufficiently  large 
to  “skip”  over  the  weak  and  rapidly  decaying  up¬ 
stream  propagation  of  information.  Alternatively,  ad¬ 
ditional  approximations  can  be  introduced  to  arrive 
at  a  well-posed  parabolic  problem  that  can  be  solved 
without  restriction  on  the  size  of  the  marching  step. 
Contrary  to  other  conclusions  [21],  well-posedness  of 
the  problem  cannot  be  enforced  by  a  proper  choice 
of  the  wavenumber  a  but  only  by  neglecting  small 
terms  to  suppress  acoustic  sources  and  even  smaller 
terms  to  suppress  viscous  sources  of  ellipticity.  The 
eflFect  of  these  steps  can  be  demonstrated  by  numer¬ 
ical  results.  Other  formulations  of  the  problem  such 
as  the  velocity-vorticity  formulation  studied  by  Haj- 
Hariri  [20]  exhibit  similar  propagation  characteristics. 
Bertolotti  (personal  communication,  1993)  has  not 
found  any  step-size  restriction  in  his  work  with  the 
Orr-Sommerfeld-Squire  formulation  [19,  22]. 


C=-eo^A^  +  (C7-|^)A-±,^+K,.(236) 

The  ellipticity  of  the  original  equations  is  completely 
removed  if  the  roots  A  =  A(/z)  of  the  determinant 
det  A  =  0  are  real  for  all  real  values  of  /z.  We  obtain 

det  A  =  (ciC2A^  -I-  /z^) 

[coA^  ^  {ReU  -  2za)A  +  /z2  ^  ReVfi]^ 

X  (24) 

and  immediately  realize  that  there  are  two  imaginary 
roots 

A  =  diz^z  ,  (25) 

except  if  either  ci  =  0,  C2  =  0,  or  both.  Paraboliza¬ 
tion,  Co  =  0,  has  no  effect  on  these  roots  of  acoustic 
origin.  To  remove  this  source  of  ellipticity,  we  can 
neglect  either  Ux  in  the  continuity  equation  or  Px  in 
the  x-momentum  equation.  The  upstream  propaga¬ 
tion  of  information  is  not  related  to  the  convective 
terms,  as  suggested  by  Chang  et  al.  [21]. 

For  the  complete  equations,  the  second  set  of  roots 
of  multiplicity  2  is 

^  =  -h{ReU-2ia) 


3.2.2  Closure  and  Norm 

Solving  eqs.  (19)  as  an  initial-boundary  value  problem 
is  not  straightforward.  The  appearance  of  dafdx  is 
disturbing  and  a  affects  the  operators,  yet  is  a  priori 
unknown.  Key  to  determining  a(x)  is  the  ambiguity 
in  the  partition  (16)  of  q  which  permits  the  incorpo¬ 
ration  of  streamwise  changes  of  the  wave  function  x 
into  q.  Given  some  Aa,  we  can  write  eq.  (16)  as  well 
as 

q'{x,y,z,t)  =  q{x,y)x{x,z,t)  ,  (28a) 

where 

q  =  qexp(zAa)  ,  x  =  xexp(— zAa)  .  (286) 

An  additional  equation  is  necessary  to  remove  this 
ambiguity.  We  note  that  Aa  affects  the  derivatives 
qx  and  q^x-  To  reduce  the  effect  of  the  PSE  approx¬ 
imation,  we  aim  at  finding  a  relation  that  removes 
the  oscillatory  behavior  from  q.  We  realize,  though, 
that  with  the  choice  of  an  appropriate  Aor,  we  can 
achieve  this  goal  only  locally  at  some  point  pm  or  in 
an  average  sense  across  the  y  domain. 

In  the  absence  of  any  information  on  streamwise 
changes,  we  can  impose  the  restriction 


otx  —  9 


(29) 
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and  absorb  all  streamwise  changes  into  q{x,y).  This 
choice  is  most  convenient  for  the  first  step  of  the  nu¬ 
merical  integration  of  eq.  (19).  However,  application 
of  eq.  (29)  for  all  x  would  maintain  the  error  ow¬ 
ing  to  the  neglect  of  higher  streamwise  derivatives  at 
0{Aa^)  and  prevent  convergence  to  the  correct  solu¬ 
tion. 

The  second  choice  imposes  a  streamwise  indepen¬ 
dent  norm  on  q.  Such  a  norm  could  be  applied  lo¬ 
cally  at  a  fixed  position  y  —  ymi^)  as  suggested  by 
Bertolotti  [19],  where  ym  is  the  position  of  the  max¬ 
imum  streamwise  rms  fluctuation  \/2|u|.  The  defor¬ 
mation  of  the  function  lu]  (especially  in  nonlinear 
problems)  may  lead  to  jumps  of  ym  and  undesirable 
effects  on  the  solution  procedure.  We,  therefore,  pre¬ 
fer  an  integral  norm  that  is  both  physically  and  math¬ 
ematically  meaningful. 

The  complex  wave  number  a  in  eq.  (10)  is  given  by 
the  logarithmic  derivative 

-i(ln  q')^:  =  ~  a  .  (30) 

q 

From  eq.  (16),  we  obtain 

-2(lnq')r  =  a  -  ,  (31) 

where  the  last  term  causes  a  dependence  of  the  local 
wave  number  on  y,  except  if  q^r  =  0.  To  remove  this 
dependency,  we  multiply  eq.  (31)  with  the  weight  |qp, 
integrate  over  the  domain  Q  in  y,  and  divide  by  the 
integral  of  |q|^.  These  steps  yield 


■  In  _  ,4_qtq^ 

fn  ^  *  in 


(32) 


where  the  dagger  denotes  the  complex  conjugate.  We 
now  choose 


i 


q^qxdy  =  0 


(33) 


to  normalize  q  and  consider  eq.  (32)  the  definition  of 
a(x). 

Eq.  (33)  minimizes  the  streamwise  change  5q/5r 
in  a  weighted  sense  across  the  y  domain.  Alternative 
conditions  can  be  formulated  using  different  weight¬ 
ing  functions  or  deleting  certain  components  of  q. 
Particularly  useful  is  the  definition  of  a  obtained  from 
eq.  (32)  by  deleting  the  pressure  and  replacing  q  by 
the  vector  v  of  velocity  amplitudes.  The  denomina¬ 
tor  then  is  the  kinetic  energy  integral,  and  the  norm 
is 

I  ^^^xdy  =  0  .  (34) 

Jci 


The  even  simpler  condition 


Ju^Uxdy  =  0 
n 


(35) 


has  been  successfully  applied  in  previous  work  [19, 

22], 


Equations  (19)  together  with  one  of  the  conditions 
(33)  to  (35),  or  any  other  suitable  norm  are  the  closed 
system  of  parabolized  stability  equations  (PSE).  This 
system  permits  the  simultaneous  calculation  of  a(x) 
and  q(a),  y)  in  a  streamwise  marching  procedure.  Dif¬ 
ferent  norms  lead  to  different  partitions  (16)  of  the  so¬ 
lution  q'  and  to  similar  yet  different  results  for  a{x) 
and  q(x*,  y).  The  physical  solution  q',  however,  is  the 
same  to  within  small  effects  of  the  norm  on  the  PSE 
approximation.  Work  is  currently  conducted  to  im¬ 
plement  an  optimal  norm  that  would  minimize  the 
effect  of  the  PSE  approximation. 

In  parallel  flows,  K  =  0,  the  choice  of  differ¬ 
ent  norms  has  no  effect  on  the  asymptotic  result 
lirur-^oo  Oi  =  a,  lima^^oo  q  =  q  since  the  wavenumber 
is  independent  of  y  and  the  same  for  all  components 
of  q.  For  parallel  flows,  it  can  be  shown  that  the  PSE 
approach  leads  to  the  correct  solution  of  the  elliptic 
problem  provided  Aa  =  a  -  oq  and  Aq  =  q  ™  qo 
are  sufficiently  small,  where  cto  and  qo  are  initial  es¬ 
timates  similar  to  those  used  for  the  iterative  solution 
of  the  local  eigenvalue  problem. 


3.2.3  Local  Analysis  of  Weakly  Nonparallel 
Flow 

All  previous  work  on  stability  theories  for  weakly  non¬ 
parallel  flows  has  reduced  the  problem  to  ordinary 
differential  equations  for  local  solutions.  These  the¬ 
ories  have  been  criticized  for  not  being  rational  per¬ 
turbation  expansions  since  the  lowest  order,  the  Orr- 
Sommerfeld  problem  for  parallel  flow,  includes  terms 
of  the  same  order  0(Re~^)  as  the  nonparallel  cor¬ 
rection.  Based  on  the  parabolized  equations,  we  can 
develop  a  local  formulation  that  does  not  exhibit  this 
shortcoming. 

Let  us  approximate  a,  q,  and  the  basic  flow  by 
their  lowest-order  Taylor  expansions  with  respect  to 
X  in  the  neighborhood  of  some  a:o: 


a(x’)  =  a(xo)  +  (x  -  xq) 


da 

dx 


laro 


:::r  rvo  4-  (x  —  Xo)^!  ,  (36a) 

q(2:>  y)  =  q(a;o,  Z/)  +  (a- -  a:o)  ^ 

^0 

=  qo(y)  +  (a;  -  a;o)qi(y)  .  (366) 

Introducing  these  expansions  into  eq.  (19)  and  requir¬ 
ing  the  equations  to  be  valid  for  varying  x  we  obtain 
two  coupled  ordinary  differential  equations, 


Lq  +  Li  -h  aiL^ 
L4  -f  aiL2 


(37) 


where  the  operators  are  obtained  from  L,M,N  in 
the  following  way:  partial  differentiation  is  replaced 
by  ordinary  differentiation;  a  is  replaced  by  oq;  Lq 
is  the  strictly  parallel  part  of  L,  Li  the  nonparallel 


4-9 


correction  of  L,  L2  —  M ^  L3  =  N ,  and  L4  originates 
from  the  expansions  of  U  and  W  in  L.  The  homo¬ 
geneous  boundary  conditions  on  the  velocities  in  q 
apply  separately  to  qo  and  qi. 

The  homogeneous  system  associated  v^^ith  eqs.  (37) 
represents  a  coupled  system  of  equations  for  the  un- 
kno^vn  quantities  ao,  ai,  qo,  and  qi  subject  to  two 
(complex)  normalizations.  The  first  norm  removes 
the  ambiguity  in  the  eigensolution.  The  second  norm 
governs  the  split  of  the  x  dependence  between  q  and 
a  as  with  eq.  (19). 

In  the  simplest  case,  we  can  apply  the  restriction 
(29),  ai  =  0,  and  solve 

subject  to  homogeneous  boundary  conditions  and  a 
norm  on  the  eigenvector.  The  modified  complex 
wavenumber  can  be  obtained  from  eq.  (31).  As  usual 
in  nonparallel  flows,  the  growth  rate  is  different  for 
different  flow  quantities  and  depends  on  the  y  loca¬ 
tion. 

As  an  alternative,  one  can  assume  that  the  solu¬ 
tion  of  the  weakly  nonparallel  problem  differs  only  by 
small  terms  from  the  solution  of  the  strictly  parallel 
problem  Loq*  =  0  and  attempt  an  iterative  solution 
of  the  system  (37).  This  procedure  requires  separate 
norms  on  qo  and  qi  and  ori  is  determined  by  a  solv¬ 
ability  condition.  The  first  iteration  loop  is  identical 
with  the  method  described  by  Saric  and  Nayfeh  [18]. 

Determining  the  solution  of  the  nonparallel  prob¬ 
lem  by  solving  eqs.  (37)  accounts  simultaneously  for 
the  terms  of  orders  C7(l)  and  0{Re'~^)  and  is,  there¬ 
fore,  not  subject  to  the  criticism  mentioned  above. 
The  direct  relation  to  the  method  of  Saric  &  Nayfeh 
shows  on  the  other  hand,  that  accounting  for  the 
terms  of  0{Re~^)  at  different  orders  of  approxima¬ 
tion  is  no  shortcoming  of  their  method. 

The  local  analysis  of  weakly  nonparallel  flows  is 
applicable  for  waves  with  sufficiently  large  wavenum¬ 
bers,  i.e.  if  the  changes  of  the  basic  flow  over  one 
wavelength  are  sufficiently  small.  The  procedure 
breaks  down  for  long  waves  as  they  appear  at  low 
Reynolds  numbers  or  low  frequencies  [22].  The  parab¬ 
olized  stability  equations  are  well  suited  to  analyze 
long  and  short  waves. 

3.2.4  Some  Formal  Properties 

Convergence  of  the  PSE  solution  to  the  correct  re¬ 
sult  is  nontrivial  and  cannot  be  expected  for  all  sta¬ 
bility  problems.  An  obvious  counter-example  is  the 
streamwise  evolution  of  an  initial  disturbance  in  cir¬ 
cular  Couette  flow,  where,  after  integrating  over  the 
circumference,  the  condition  of  periodicity  would  be 
violated.  The  neglect  of  the  second  and  higher  deriva¬ 
tives  of  a  and  q  with  respect  to  x  and  the  change  of 
the  mathematical  character  of  the  governing  partial 


differential  equations  from  elliptic  to  quasi-parabolic 
is  only  permitted  if  the  stability  problem  is  governed 
by  downstream  propagating  information  while  the 
upstream  propagation  can  be  neglected. 

The  use  of  parabolic  differential  equations  for  ana¬ 
lyzing  problems  of  basically  elliptic  nature  with  small 
feedback  is  successful  in  some  other  areas,  e.g.  in  the 
analysis  of  acoustic  wave  propagation  [23].  In  the 
field  of  flow  instabilities,  the  propagation  character¬ 
istics  of  disturbances  have  been  studied  in  detail  and 
led  to  the  classification  into  absolute  and  convective 
instabilities  [12].  In  these  terms,  the  PSE  approach  is 
valid  for  convectively  unstable  flows.  This  class  con¬ 
tains  pipe  and  channel  flows  and  the  technologically 
important  boundary  layers,  mixing  layers,  far  wakes, 
and  others  that  make  it  worthwhile  pursuing  this  new 
avenue. 

The  convective  nature  of  linear  instabilities  can  be 
assured  by  analytical  studies.  For  nonlinear  prob¬ 
lems  such  as  the  later  stages  of  transition,  the  role  of 
upstream  propagating  information  can  only  be  esti¬ 
mated  from  experimental  or  numerical  evidence.  Our 
applications  of  the  PSE  approach  suggest  that  the 
transition  process  in  boundary  layers  is  largely  gov¬ 
erned  by  downstream  propagating  information. 

In  contrast  to  the  normal- mode  approach,  the  PSE 
formulation  and  marching  solution  neither  exploits 
nor  requires  the  homogeneity  of  the  problem.  This 
property  permits  both  inhomogeneous  boundary  con¬ 
ditions  and  a  forcing  function  instead  of  zero  on 
the  right  hand  side  of  eq.  (19).  The  approach  can, 
therefore,  handle  forced  and  natural  oscillations  of 
the  flow  within  a  largely  unified  numerical  approach. 
Receptivity,  the  reaction  of  the  basic  flow  to  small 
forced  disturbances,  is  a  yet  insufficiently  explored 
field  of  transition  research  which  can  benefit  from 
the  PSE  approach.  Certain  types  of  bypass  transi¬ 
tion  that  are  associated  with  strong  transient  energy 
growth  although  linear  stability  would  predict  stabil¬ 
ity  [24,  25,  26,  27]  are  directly  accessible  to  the  PSE 
approach.  The  forcing  function  can  also  describe  the 
nonlinear  coupling  in  the  simultaneous  evolution  of  a 
group  of  modes  distinguished  by  their  frequencies  and 
spanwise  wavenumbers.  These  properties  of  the  PSE 
permit  the  analysis  of  transition  from  the  ingestion 
of  small  disturbances  through  primary  and  higher  in¬ 
stabilities  to  the  breakdown  of  the  laminar  basic  flow. 
Together  with  efficient  marching  techniques,  the  PSE 
appear  as  a  suitable  basis  for  applications  in  both  re¬ 
search  and  engineering  practice. 

3.2.5  Solution  Procedure 

Before  solving  eqs.  (19)  and  (34),  we  need  to  spec¬ 
ify  initial  and  boundary  conditions.  For  comparison 
with  the  local  analysis  in  section  3.1.1,  we  specify  Re, 
(3,  cj,  and  V,  provide  the  same  estimates  cvq  for  the 
eigenvalue  and  qo  for  the  eigenfunction  as  initial  con¬ 
ditions  at  x’5,  and  apply  the  homogeneous  boundary 


Lq  4-  L\  L2 
L4  Lq 
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conditions  (11). 

To  solve  the  PSE,  the  diflerential  operators  (21) 
and  boundary  conditions  are  converted  into  algebraic 
form  by  any  method  suitable  for  the  Orr-Sommerfeld 
problem.  We  use  a  spectral  collocation  method 
with  Chebyshev  polynomials,  together  with  a  domain 
transformation  if  appropriate.  Alternatively,  we  use 
two-point  or  three-point  Hermitian  high-order  finite- 
difference  methods  on  nonuniform  grids.  These  meth¬ 
ods  lead  to  algebraic  systems  with  banded  or  block- 
tridiagonal  matrices  which  can  be  more  efficiently 
solved  for  the  relatively  large  system  (18)  in  natural 
variables. 

The  streamwise  derivative  is  in  the  simplest  case 
approximated  by  the  two-point  finite  difference 

<ir «  -qj) .  (39) 

where  J  >  0  is  the  step  index  and  q^-  ^j{y)  — 

q(xj ,  y).  Since  a  appears  in  the  differential  operators, 
the  system  is  nonlinear,  and  a  predictor-corrector  ap¬ 
proach  is  employed  to  find  the  solution  at  Xj^i.  The 
explicit  forward  difference  formula  leads  to  numeri¬ 
cal  instability.  Although  second-order  accurate,  the 
central  difference  formula  which  evaluates  qa;  at  the 
midpoint  Xj^if2  has  not  been  found  advantageous 
in  comparison  with  the  first-order  accurate  backward 
difference 

(qr)j+i «  ^(qj+i  -  qj)  •  (40) 

Equation  (19)  then  assumes  the  form 


=  ^iViq;-  .  (41) 

where  (ar)j4-i  =  ““  l^he  index  n  counts 

the  iterations.  Starting  with  ==  ckj  according 

to  eq.  (29)  or  with  any  better  initial  estimate,  we 
obtain  q!j^i  and  henceforth  exploit  eq.  (32)  to  obtain 
an  updated  wavenumber 


^  n+l  n 

—  ^j  +  l 


.  (42) 


The  integration-update  cycle  is  repeated  until  eq. 
(33)  is  satisfied  within  a  given  error  limit  and  the 
solution  has  converged  to  ay+i,  qj_}.i(y).  We  then 
proceed  to  the  next  step  in  x  where  a  better  estimate 
for  a  and  is  known.  The  marching  procedure  ter¬ 
minates  if  the  convergence  criterion  is  not  satisfied 
within  a  given  number  of  iterations. 

A  simultaneous  iteration  of  q”^^  and  can  be 
implemented  e.g.  by  using  Newton’s  method  on  the 
coupled  equations  (41)  and  (42). 


3.2.6  Scaling  Considerations 


Similar  to  DNS,  the  solution  of  the  PSE  refers  to  fixed 
scales  Uq  and  =  <^*(0:0)  for  velocity  and  length, 
where  6;{x)  ^  and  a:  =  x*/6*  (the 

star  denotes  dimensional  quantities).  The  Reynolds 
number  Reo  =  is  fixed.  The  reference  po¬ 

sition  xq  is  not  necessarily  the  initial  position 
In  contrcLst,  local  stability  analyses  commonly  use 
the  local  length  scale  S*  and  the  Reynolds  numbers 

R^x  =  Re  ~  For  example,  the 

frequency  a;  =  in  the  PSE  formulation  de¬ 

scribes  a  wave  of  fixed  physical  frequency,  and  there 
is  no  need  to  refer  to  the  usual  nondimensional  fre¬ 
quency  F.  To  avoid  confusion  between  the  differently 
scaled  quantities,  we  use  the  caret  to  indicate  local 
scales,  e.g.  cj  =  u}*6*{x)/U*. 

For  Blasius  flow,  the  two  scalings  provide  F  = 
uj/Reo  =  u}/Re  and  the  relations 


Re  _  Re  ^  Re  ^  ^ 

Since  x  is  scaled  with  6q,  we  also  have  the  relations 


Rej;  =  xReo  ,  Re  -  y/xR^  .  (44) 

Similar  relations  can  be  established  for  other  basic 
flows. 

In  nonparallel  flow,  additional  corrections  to  a 
may  be  necessary  to  account  for  specific  quantities 
and/or  positions  y.  The  results  reported  in  the  fol¬ 
lowing  are  based  on  the  norm  (34)  and  the  ampli¬ 
tude  A  =  (f^  |v[^d2/)^/^.  For  comparison  with  exper¬ 
iments,  we  convert  to  the  local  magnitude  of  velocity 
components  (for  steady  modes)  or  the  rms  value  (for 
oscillatory  modes). 

The  nature  of  the  primary  results  obtained  by  the 
PSE  approach  and  linear  stability  theory  (LST)  is 
different  and  their  comparison  involves  proper  con¬ 
version.  Since  the  PSE  results  for  the  flow  field  in 
space  and  time  (similar  to  DNS  solutions)  are  closer 
to  the  answers  sought  by  the  analysis  than  eigenval¬ 
ues  as  a  function  of  the  local  Reynolds  number,  we 
prefer  the  fixed  scaling  and  present  the  evolution  over 
X  or  Rcx  according  to  eq.  (44). 


3.3  Some  Linear  Studies 

As  an  introduction  to  PSE  applications,  we  consider 
some  applications  of  the  linear  equations.  Thes  ap¬ 
plications  concern  both  instability  and  receptivity  of 
2D  and  3D  boundary  layers. 

As  a  generic  example,  we  first  analyze  the  2D 
boundary-layer  flow  over  a  flat  plate  at  zero  pressure 
gradient.  The  examples  can  be  considered  as  a  “test 
suite”  for  the  development  of  a  PSE  code. 


3,3.1  Strictly  Parallel  Blasius  Flow 

As  a  first  test  of  the  PSE  concept,  we  apply  the  solu¬ 
tion  procedure  to  a  strictly  parallel  flow  which  is  in- 


4-11 


dependent  of  x.  We  choose  the  Blasius  profile  and  a 
two-dimensional  mode  with  lj  =  0.0344  at  Rcq  =  400 
(F  =  86  •  10^®).  The  basic  flow  is  held  fixed,  indepen¬ 
dent  of  X.  We  start  the  computation  at  ~  xq  =  400 
and  choose  as  an  initial  condition  the  results  of  the 
local  analysis  (a)  for  the  same  conditions,  (b)  for 
u)  =  0.0172  at  the  same  Reynolds  number.  Giving  the 
value  of  a  for  case  (b)  as  an  initial  guess  for  the  local 
analysis  at  a;  =  0.0344  leads  to  convergence  to  some 
other  eigenvalue  with  strong  damping.  Indices  s  and 
0  distinguish  the  starting  position  for  the  PSE  run 
and  the  reference  position,  respectively.  The  results 
of  the  PSE  calculation  for  the  streamwise  evolution 
of  ar  and  ai  are  shown  in  figure  2. 


Figure  2:  Evolution  of  the  wave  number  a  for  differ¬ 
ent  initial  conditions. 

Case  (a)  (dashed)  is  a  valid  test  of  the  integrity  of 
the  PSE  code  which  should  maintain  the  correct  re¬ 
sult  over  the  marching  domain.  The  other  case  serves 
to  demonstrate  the  proper  function  of  the  wavenum¬ 
ber  update  by  using  a  suitable  norm.  For  initial 
condition  (b),  the  first  few  marching  steps  provide  a 
rapidly  decaying  transient  and  the  PSE  result  asymp¬ 
totes  to  exactly  the  same  value  of  a  obtained  with 
the  local  analysis  at  comparable  numerical  resolu¬ 
tion.  As  an  advantage  we  recognize  that  even  poor 
initial  conditions  lead  to  convergence  to  the  most  un¬ 
stable  mode.  Since  \ai\  is  typically  much  smaller 
than  ar,  the  initial  transients  of  the  growth  rate  ap¬ 
pear  stronger.  The  results  are  largely  insensitive  to 
changes  of  the  marching  step  Ax  and  the  choice  of  al¬ 
ternative  integration  schemes.  The  transient  solution 
closely  approximates  the  physical  behavior  provided 
the  initial  equation  error  is  sufficiently  small. 

Since  neither  the  derivation  nor  the  solution  pro¬ 
cedure  for  the  PSE  requires  the  problem  to  be 


homogeneous,  we  can  introduce  disturbances  not 
only  through  the  initial  conditions  but  also  through 
boundary  conditions  or  forcing  inside  the  flow  do¬ 
main,  as  in  laboratory  experiments.  In  the  second 
test,  we  leave  the  initial  conditions  identically  zero 
and  apply  a  small  v  component  at  the  wall  to  gen¬ 
erate  a  disturbance  of  frequency  u)  =  0.0344  at  the 
same  values  of  i?eo  and  x,  as  before,  while  =  0 
for  X  >  x^.  Figure  3  shows  the  result  for  the  evolu¬ 
tion  of  a  in  this  case.  After  a  short  period  of  “wild” 
transient  behavior,  the  solution  settles  at  the  same 
values  given  by  the  previous  runs  and  the  local  anal¬ 
ysis.  The  capability  to  force  disturbances  in  the  ab¬ 
sence  or  presence  of  initial  conditions  is  key  to  many 
studies  on  receptivity  and  transition. 


Figure  3:  Disturbance  introduced  by  a  small  v  com¬ 
ponent  at  the  initial  position  x,. 

The  application  of  the  PSE  to  parallel  flows  is  a 
demonstration  of  the  concept  and  a  verification  of  the 
stability  and  robustness  of  the  marching  procedure. 
Considering  the  effort  for  marching  with  a  step  size 
of  Ax  «  20  over  the  range  of  x  shown  in  figures  2 
and  3,  the  direct  solution  of  the  eigenvalue  problem 
is  indeed  more  efficient.  However,  this  comparison 
changes  drastically  if  the  basic  flow  is  nonparallel  and 
the  coefficients  of  the  stability  equations  change  with 
the  marching  variable  x. 

3.3.2  Nonparallel  Blasius  Flow 

To  illustrate  some  advantages  of  the  PSE  technique, 
we  consider  the  evolution  of  a  two-dimensional  wave 
at  frequency  uj  =  0.0344  with  Reo  =  400  for  x  > 
X,  =  400.  The  initial  conditions  are  given  by  the 
local  analysis  at  x,. 

Figure  4  shows  the  wavenumbers  and  growth  rates 
obtained  from  the  local  analysis  in  comparison  with 


Figure  4:  Evolution  of  the  wavenumber  a  in  a  Blasius 
boundary  layer.  LST  (dashed)  and  PSE  (full  line). 

the  result  of  the  PSE.  With  LST  we  refer  to  the  tra¬ 
ditional  linear  stability  theory  for  locally  parallel  flow 
while  the  “modified  LST”  includes  the  transverse  ve¬ 
locity  V  and  the  derivative  Ux- 

We  observe  only  small  differences  in  the  growth 
rates.  We  also  find  that  the  variation  of  a  is  very 
small,  as  reported  by  Mack  [13],  and  the  effect  of 
is  negligible.  The  variation  of  d  =  aRe/Reo  in  results 
of  the  traditional  analysis  is  largely  due  to  the  stream- 
wise  variation  of  the  length  scale  (5;!:(ar)  while  the  PSE 
results  refer  to  the  fixed  length  scale  6q  =  (5*(a:o). 
The  jiggles  near  the  starting  position  are  caused  by 
the  marginal  step  size  of  Ax  =  10  and  disappear  as 
decreases. 

To  illustrate  the  changes  in  the  amplitude  functions 
over  the  large  x  domain  in  figure  4,  figure  5  compares 
the  normalized  streamwise  rms  fluctuations 
at  different  positions.  The  modest  changes  justify  the 
neglect  of  the  derivatives  in  the  PSE. 

For  the  numerical  work,  it  is  important  to  note 
that  the  boundary  layer  actually  grows  in  the  Carte¬ 
sian  coordinates  used  here.  The  numerical  grid  and 
its  extent  to  some  ymax  must  be  carefully  chosen  to 
cover  the  boundary  layer  at  the  end  of  the  run  yet  to 
maintain  sufficient  resolution  for  the  amplitude  func¬ 
tions  of  the  disturbance  near  the  wall. 

In  general,  the  results  of  the  PSE  runs  are  very 
sensitive  to  small  changes  in  basic  flow,  initial  con¬ 
ditions,  and  numerical  treatment.  This  property  of¬ 
ten  makes  it  difficult  to  find  agreement  between  re¬ 
sults  of  different  codes  or  different  authors  for  vir¬ 
tually  the  same  case.  This  fact  is  especially  true 
for  the  amplitude-growth  curves  which  are  more  rel¬ 
evant  in  practical  applications  than  the  growth  rates 
provided  by  the  local  analysis.  Small  deviations  in 


Figure  5:  Normalized  rms  profiles  of  the  amplitude 
function  u  at  x  =  500  (dashed)  and  2000  (full  line). 

growth  rates  accumulate  by  integration  over  an  ex¬ 
tended  streamwise  distance.  The  PSE  approach  pro¬ 
vides  these  amplitude-growth  curves  directly  together 
with  ol  and  the  spatial  growth  rates 

The  previous  comparison  of  local  and  PSE  results 
for  two-dimensional  waves  seems  to  indicate  a  rather 
small  effect  of  the  nonparallelism.  This  conclusion 
changes,  however,  for  oblique  or  three-dimensional 
waves  with  {5 


Figure  6:  Comparison  of  results  for  the  wavenumber 
of  a  three-dimensional  disturbance.  LST  (dashed) 
and  PSE  (full  line). 

The  streamwise  variation  of  a  for  a  wave  with 
Lj  =  0.02604  and  j3  =  0.14  at  Re^  =  400  (the  sub¬ 
harmonic  wave  in  the  Kachanov-Levchenko  experi¬ 
ment  [10])  is  shown  in  figure  6.  Starting  from  the 
initial  value  given  by  the  local  analysis,  the  growth 
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rate  develops  to  a  more  unstable  level  which  it  main¬ 
tains  throughout  the  run.  For  smaller  j3,  the  PSE 
analysis  exhibits  instability  while  the  wave  remains 
stable  according  to  the  locally  parallel  theory. 

The  local  theory  for  weakly  nonparallel  flows  of  sec¬ 
tion  3.2.3  can  retrieve  some  of  the  effects  of  boundary- 
layer  growth,  provided  the  wavenumber  is  not  too 
small.  As  the  wave  angle  increases,  this  condition  is 
not  satisfied  and  the  local  theory  fails.  This  failure 
has  been  shown  by  Hall  [28]  for  the  case  of  Gortier 
vortices  (wave  angle  90®).  Hall  demonstrated  that 
Gortier  vortices  are  governed  by  an  initial  value  prob¬ 
lem  for  parabolic  stability  equations.  These  equa¬ 
tions  are  derived  from  the  boundary-layer  equations, 
and  therefore,  the  pressure  terms  p  and  px  in  the  x- 
momentum  equation  (18b)  are  absent.  Since  Gortier 
vortices  are  steady,  the  analysis  can  set  or  =:  0  and 
absorb  the  growth  into  the  amplitude  function. 

3.3.3  Vortices  in  Blasius  Flow 

In  the  flat-plate  boundary  layer,  steady  disturbances 
have  not  been  revealed  by  the  local  stability  anal¬ 
ysis  although  quasi-steady  regular  spanwise  modu¬ 
lations  of  the  basic  flow  have  been  observed  in  nu¬ 
merous  experiments  since  Dryden  [29]  and  studied  in 
detail  [30,  31,  32].  Even  in  wind  tunnels  with  low 
turbulence  levels,  these  “Klebanoff  modes”  can  grow 
with  «  to  amplitudes  of  5%  before  the  onset 
of  TS  instability.  Similar  modulations  have  been  ob¬ 
served  in  channel  flow  [33]  where  the  amplitude  varies 
slowly.  The  modulations  in  channel  flow  are  associ¬ 
ated  with  the  stable  yet  very  weakly  damped  longi¬ 
tudinal  vortices  and  tz-modes  [34]  that  are  eigenso- 
lutions  of  the  Orr-Sommerfeld  and  Squire  equation, 
respectively,  for  a  =  0,  0.  In  the  Blasius  flow, 

these  modes  are  hidden  in  the  continuous  spectrum 
and  the  local  analysis  cannot  provide  initial  data  for 
the  marching  technique. 

With  proper  boundary  conditions,  the  PSE  analy¬ 
sis  for  steady  modes,  a;  =  0,  can  be  performed  with 
zero  initial  conditions.  To  introduce  u-modes,  which 
in  a  parallel  flow  have  t;  =  ly  =  0,  disturbances  can 
be  generated  by  a  small,  spanwise  periodic  variation 
of  the  wall-shear  stress,  Du  =  e  at  y  =  0,  ar  = 

Test  runs  show  the  rapid  formation  of  a  streamwise 
u  and  weaker  v  and  w  components  that  decay  with 
cti  =  0(10“'^).  Owing  to  the  nonparallelism  of  the 
basic  flow  and  the  nonzero  derivative  Ua?,  pure  u- 
modes  cannot  appear  in  the  Blasius  flow.  Since  nei¬ 
ther  the  streamwise  amplitude  variation  nor  the  u 
profile  is  consistent  with  observations,  the  spanwise 
modulations  are  likely  associated  with  other  mecha¬ 
nisms. 

To  study  the  generation  and  evolution  of  forced 
longitudinal  vortices,  disturbances  can  be  generated 
by  a  small,  spanwise  periodic  variation  of  the  nor¬ 
mal  velocity  component  at  the  wall,  v(0)  =  6  at 
X  =  Xg.  Figure  7  shows  the  evolution  of  the  am- 
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Figure  7:  Growth  of  longitudinal  vortices  forced  at 
different  initial  positions. 

plitude  for  three  initial  positions,  x^  =  25,100,400, 
with  P  =  0.14,  RtQ  =  400,  and  e  =  10“^.  The  span- 
wise  wavenumber  is  larger  than  the  value  (5  »  0.105 
(/?  0.25  at  Re  =  941)  observed  by  Klebanoff  et 

al.  [8]  and  about  half  the  value  /?  0.33  ( /?  0.6  at 

Re  =  717)  reported  by  Kendall  [35,  32].  The  small 
initial  disturbances  monotonically  grow  to  consider¬ 
able  amplitudes.  The  growth  rate  increases  with  x,. 
The  velocity  profiles  of  the  disturbances  closely  re¬ 
semble  those  of  Gortier  vortices,  with  a  pronounced 
u  component  and  v  and  w  typically  two  orders  of 
magnitude  smaller.  The  amplitude  evolution  of  dis¬ 
turbances  with  different  spanwise  wave  numbers  is 
given  in  figure  8  for  x,  =  25  and  the  previous  values 
of  i2eo  and  e.  The  curves  for  the  larger  wavenum¬ 
bers  show  that  the  initial  growth  ultimately  changes 
into  decay.  The  maximum  amplitude  is  higher  and 
appears  at  lower  x  as  f3  increases. 


Figure  8:  Growth  curves  of  forced  longitudinal  vor¬ 
tices  for  different  spanwise  wavenumbers 

At  X  =  1285,  Re  =  717,  the  position  of  Kendall’s 
detailed  measurements,  the  largest  amplitude  occurs 
for  /?  «  0.5,  in  the  same  range  as  in  the  experiments. 
The  profile  of  the  calculated  u  disturbance  shown  in 
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figure  9  with  its  maximum  near  r]  -  yReo/ Re  ~  2.2 
is  similar  to  the  observed  distribution  (crosses)  in  fig¬ 
ure  10.  The  accompanying  curve  would  result  from 
a  spanwise  thickening  and  thinning  of  the  bound¬ 
ary  layer,  as  pointed  out  by  Taylor  [36]  and  Kle- 
banoff  [30]. 


Figure  9:  Calculated  velocity  distribution  of  the  Kle- 
banoff  mode. 


FLUCTUATION  PROFILES  WITH  GRID  TURBULENCE 


Figure  10:  Meetsured  velocity  distribution  of  the  Kle- 
banoff  mode. 

The  single  disturbance  introduced  at  Re  —  25  de¬ 
cays  at  the  position  Re  =  717  of  the  observation.  In 
the  experiment,  however,  the  disturbances  are  con¬ 
tinuously  seeded  by  area-distributed  receptivity  to  a 
wide  spectrum  of  wavenumbers.  To  better  reproduce 
the  experimental  conditions,  a  discrete  spectrum  of 
modes  in  z  and  proper  inhomogeneous  boundary  con¬ 
ditions  in  some  range  of  x  should  be  applied.  Be¬ 
sides  through  i;(0),  the  Klebanoff  modes  can  be  intro¬ 
duced  through  inhomogeneous  conditions  on  Dw  at 
the  wall,  and  most  efficiently,  through  spanwise  vari¬ 
ations  of  the  wall  pressure  p(0).  While  the  strongest 
response  is  for  steady  modes,  the  mechanism  covers 
a  broad  band  of  low  frequencies.  For  a  smooth  plate, 


generation  of  Klebanoff  modes  by  spanwise  pressure 
variations  associated  with  free-stream  turbulence  is 
the  most  likely  scenario. 

The  significant  transient  amplitude  (or  energy) 
growth  of  steady  and  low-frequency  disturbances  in 
the  absence  of  linear  instability  is  an  example  of  the 
bypass  mechanism  described  by  Gustavsson  [24]  and 
Henningson  [25]  for  channel  flow.  Related  studies  of 
optimal  energy  growth  in  a  crude  model  of  Blasius 
flow  by  Butler  &  Farrel  [26]  show  that  the  optimum 
is  associated  with  steady  disturbances  of  a  spanwise 
wavelength  similar  to  those  found  in  the  experiments 
of  Kendall  and  in  our  PSE  analysis  of  forced  dis¬ 
turbances.  While  previous  studies  of  the  transient 
growth  mechanism  were  restricted  to  the  temporal 
evolution,  the  PSE  are  well  suited  to  analyze  the  spa¬ 
tial  evolution  and  properly  account  for  the  nonparal¬ 
lelism  of  the  flow. 

3.3.4  Effect  of  Curvature 

The  similarity  between  the  velocity  distributions  of 
Klebanoff  modes  and  Gdrtler  vortices  suggest  that 
they  may  be  strongly  affected  by  curvature.  Also,  the 
occurrence  of  Gortler  vortices  on  concave  walls  may 
be  related  to  the  same  receptivity  mechanisms.  To 
analyze  the  effect  of  curvature,  the  PSE  have  been  de¬ 
rived  from  the  Navier-Stokes  equations  for  a  surface- 
oriented  coordinate  system  with  streamwise  varying 
radius  of  curvature  R{x)  [37].  While  these  equations 
permit  interesting  studies  on  the  disturbance  evolu¬ 
tion  on  wavy  walls  or  airfoils,  we  consider  here  walls 
with  a  constant  (dimensional)  radius  of  curvature  and 
K  =  Sq/R{x)  =  const,  with  respect  to  the  reference 
length  at  Rcq  —  400.  As  before,  vye  introduce  distur¬ 
bances  through  a  small,  spanwise  periodic  wall  veloc¬ 
ity  i;(0)  =  10“^  at  —  25. 


Figure  11:  Growth  curves  of  forced  longitudinal  vor¬ 
tices  for  different  convex  and  concave  wall  curvature 
in  comparison  with  the  result  for  a  flat  wall. 

Figure  11  shows  the  growth  curves  for  values  of 
K  •  10®  =  —25(10)25  together  with  the  result  for  the 
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flat  plate,  «  =  0,  for  /?  =  0.36.  After  initial  growth, 
the  amplitudes  for  /c  >  — 5  •  10“^  decrease  with  in¬ 
creasing  X,  the  amplitude  for  /c  =  — 5  •  10“®  remains 
nearly  constant,  while  the  amplitudes  for  larger  con¬ 
cave  curvature  (k  <  0)  increase  with  varying  rate. 


Figure  12:  Growth  rate  of  forced  longitudinal  vortices 
vs.  curvature  at  x  =:  1600. 

The  growth  rates  as  a  function  of  k  are  given  in 
figure  12  for  x  =  1600.  The  data  from  PSE  runs  are 
shown  by  the  circles.  The  growth  rate  varies  mono- 
tonically  as  k  increases.  After  rescaling  the  data  for 
K  =  —25  •  10^®  to  the  local  conditions  at  Re  =  800, 
the  growth  rate  d,-  =  -0.002386  agrees  surprisingly 
well  with  the  value  d,-  =  —0.002382  of  the  local  anal¬ 
ysis  of  a  Gortler  vortex  at  the  proper  parameters 
P  =  0.72,  /c  =  — 5  •  10”^.  Some  of  this  agreement 
is  incidental  since  a  PSE  analysis  for  locally  parallel 
flow,  which  is  more  consistent  with  the  local  analysis, 
provides  d,-  =  —0.002056.  However,  growth  charac¬ 
teristics  and  comparison  of  velocity  profiles  leave  no 
doubt  that  the  forced  vortices  at  «  =  —25  •  10““®  are 
Gortler  vortices.  The  PSE  analysis  clearly  shows  the 
existence  of  these  vortices  in  the  form  of  Klebanoff 
modes  on  a  flat  plate  and  increasingly  damped  modes 
on  convex  walls  of  increasing  curvature.  The  stabiliz¬ 
ing  effect  of  convex  curvature  on  longitudinal  vortices 
is  a  new  result  outside  the  scope  of  the  linear  stabil¬ 
ity  theory.  In  spite  of  the  close  relation  between  Kle¬ 
banoff  modes  and  Gortler  vortices,  it  is  appropriate 
to  maintain  the  different  names  since  Gortler  vortices 
are  inherently  associated  with  centrifugal  instability. 

3.3.5  Wave  Packets 

Besides  the  single  modes  studied  above,  the  PSE  code 
allows  tracking  the  evolution  of  waves  distinguished 
by  their  frequency  and  spanwise  wave  number.  In 
a  linear  framework,  the  waves  can  be  analyzed  se¬ 
quentially  or  simultaneously.  The  latter  approach  has 
been  chosen  since  it  enables  easier  post-processing 
of  the  data.  Test  runs  with  coarsely  spaced  modes 
have  been  performed  for  two-dimensional  wave  pack¬ 
ets  [38],  a  periodic  point  source  [39],  and  a  three¬ 


dimensional  wave  packet  [40].  More  detailed  stud¬ 
ies  with  higher  resolution  are  performed  in  cooper¬ 
ation  with  L.  Mack  to  compare  with  experimental 
and  earlier  theoretical  results.  Figure  13  shows  the 
wave  pattern  created  by  a  harmonic  point  source  at 
a:,  =  xo  =  475  with  w  =  0.0285  (i^  =  60  •  10'®). 


Figure  13:  Contours  of  equal  streamwise  disturbance 
velocity  for  0  <  x  —  x^  <  1000  at  y  =  2.5. 

The  initial  disturbance  is  modeled  by  201  modes 
with  equal  initial  amplitudes  and  phases  and  hm  — 
0.0075  •  m,  m  =  0,  •  *  • ,  200.  Alternatively,  the  initial 
conditions  could  be  derived  from  a  Fourier  decompo¬ 
sition  of  the  disturbance  source  at  the  initial  position 
by  applying  proper  boundary  conditions  and  includ¬ 
ing  selective  receptivity  mechanisms.  The  PSE  anal¬ 
ysis  is  time  consuming  yet  otherwise  straightforward 
since  the  physical  constraints  during  the  downstream 
evolution  of  the  waves  are  maintained  by  the  gov¬ 
erning  equations.  The  analysis  does  not  require  the 
heuristic  corrections  for  the  disturbance-energy  dis¬ 
tribution  across  the  boundary  layer.  Accounting  for 
a  larger  number  of  modes  is  only  a  matter  of  patience 
and  disk  space. 

3.3.6  Three-Dimensional  Boundary  Layers 

The  examples  so  far  were  for  two-dimensional  bound¬ 
ary  layers.  The  parabolized  equations  (19)  can  be  di¬ 
rectly  applied  to  three-dimensional  boundary  layers 
in  Cartesian  coordinates,  such  as  the  Falkner-Skan- 
Cooke  flows  over  a  swept  wedge  or  the  swept  Hiemenz 
flow  along  and  towards  the  sides  of  a  stagnation  line. 
For  other  cases  such  as  the  boundary  layer  over  a 
rotating  disk,  the  PSE  can  be  reformulated  in  cylin¬ 
drical  coordinates  [41].  Alternatively,  a  coordinate 
transformation  between  Cartesian  and  cylindrical  co¬ 
ordinates  can  be  prescribed  and  the  proper  equations 
can  be  generated  during  the  computation.  This  lat¬ 
ter  procedure  permits  a  modular  split  of  the  metric 
terms  from  the  formulation  of  the  PSE. 

For  the  flow  over  a  rotating  disk,  the  nondimen- 
sional  radius  is  used  as  the  marching  variable.  Fig¬ 
ure  14  shows  the  variation  of  the  growth  rate  for  a 
disturbance  with  =  30,  where  is  the  number  of 
vortices  over  the  circumference.  The  result  is  similar 
to  the  data  presented  by  Malik  k  Balakumar  [42,  fig. 
la]  who  also  compare  with  results  of  various  other 
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Figure  14:  Wave  number  and  growth  rate  of  cross- 
flow  vortices  with  np  =  30  in  the  boundary  layer  on 
a  rotating  disk. 

approaches  to  the  controversial  stability  of  this  flow. 
Additional  discussion  of  the  results  can  be  found  in 
ref.  [43].  The  figure  also  shows  some  oscillations  at 
high  Reynolds  number  that  ultimately  terminate  the 
run.  Although  the  wave  number  ar  decreases,  the  os¬ 
cillations  are  caused  by  a  strong  numerical  instability 
for  the  step  size  used  for  this  run.  With  larger  step 
size  or  control  of  the  elliptic  influence,  the  run  can  be 
continued  to  higher  Reynolds  numbers. 

The  originally  strong  interest  in  the  stability  of  the 
boundary  layer  over  a  rotating  disk  as  a  prototype 
of  boundary  layers  on  realistic  swept  wings  has  faded 
since  the  stability  properties  are  in  fact  more  com¬ 
plex  than  in  swept-wing  flows.  Falkner-Skan-Cooke 
and  swept  Hiemenz  flow  are  more  closely  related  to 
different  aspects  of  swept-wing  flows  and  we  present 
only  a  few  results  for  the  latter  case.  Based  on  eqs. 
(19),  the  marching  variable  is  x.  There  is  no  need  to 
work  in  coordinate  systems  that  adapt  one  axis  to  the 
inviscid  streamline  direction,  as  it  is  common  in  tradi¬ 
tional  N  factor  calculations.  There  is  also  no  problem 
with  the  controversial  in-plane  curvature  terms  which 
have  been  the  topic  of  various  papers  [43,  44].  Given 
the  wavenumber  jS  in  the  spanwise  z  direction,  the 
streamwise  wavenumber  ocr  adapts  properly  to  the 
local  orientation  of  the  disturbances.  Since  the  basic 
flow  is  independent  of  2r,  /?  must  be  constant  during 
the  run  [14]. 

Although  Hiemenz  flow  has  a  constant  boundary- 
layer  thickness,  the  flow  is  non-parallel,  as  shown 
by  the  streamlines.  Figure  15  shows  the  amplitude 
growth  of  steady  cross-flow  vortices  in  the  swept 
Hiemenz  flow  with  Re  =  500,  =  -0.4,  a  case  stud¬ 


Figure  15:  Amplitude  growth  curves  {N  factors)  for 
swept  Hiemenz  flow.  PSE  results  with  and  without 
(dashed)  accounting  for  nonparallel  terms. 

ied  with  DNS  by  Spalart  [45].  Re  is  formed  with 
the  spanwise  velocity  Wg  and  the  reference  length 
where  [/e(x)  —  Ax.  The  marching  vari¬ 
able  X  is  the  Reynolds  number  formed  with  f/g.  In 
general,  the  linear  stability  characteristics  provided 
by  the  local  theory  are  in  good  agreement  with  the 
PSE  results.  Accounting  for  the  variation  of  the  basic 
state  causes  some  differences  but  no  dramatic  change. 
We  note,  however,  that  these  observations  cannot  be 
generalized  to  other,  more  realistic  flows. 


Figure  16:  Wavenumber  and  growth  rate  of  steady 
cross-flow  vortices  in  swept  Hiemenz  flow  according 
to  LST  (dashed)  and  PSE  approach. 

The  wave  numbers  ar  and  growth  rates  a,-  obtained 
by  linear  theory  and  PSE  analysis  are  compared  in 


figure  16.  Similar  results  have  been  presented  by  Ma¬ 
lik  Li  [46,  47]. 


3.3.7  Receptivity  of  Hiemenz  Flow 

While  the  PSE  approach  allows  a  more  efficient,  ac¬ 
curate,  and  reliable  analysis  of  growth  rates  and  N 
factors,  it  also  provides  new  capabilities  for  ana¬ 
lyzing  and  understanding  the  receptivity  of  three- 
dimensional  boundary  layers.  Fedorov  [48]  stud¬ 
ied  the  excitation  of  cross-flow  vortices  in  subsonic 
boundary  layers  over  swept  wings  and  found  that  mi¬ 
croroughness  of  the  height  of  1%  of  the  displacement 
thickness  or  spanwise  periodic  suction/injection  with 
t;(0)  =  can  excite  cross-flow  vortices  with  an 

initial  amplitude  of  0.1%.  Radeztsky  et  al.  [49,  50] 
report  a  strong  effect  of  small  roughness  elements  on 
transition  in  swept-wing  flows.  Crouch[51]  extended 
his  local  nonlinear  approach  to  the  receptivity  of 
the  Falkner-Skan-Cooke  flow  to  steady  and  unsteady 
cross-flow  vortices  and  obtained  results  consistent 
with  the  experiments  of  Muller  &  Bippes[52].  There 
are  at  least  some  aspects  of  these  studies  that  can  be 
qualitatively  reproduced  by  modeling  the  roughness 
elements  by  a  disturbance  of  the  v  component  of  the 
velocity  normal  to  the  wall. 

Figure  17  shows  the  amplitude  of  forced  cross- 
flow  vortices  of  /?  =  —0.4  at  x  =  500  depending  on 
the  position  x,  of  the  v  disturbance.  In  all  cases, 
i;(0)  =  10""^  at  X,.  The  response  is  strongest  if  the 
disturbance  (or  the  roughness  in  the  experiments)  is 
placed  shortly  upstream  of  the  neutral  point  for  the 
onset  of  instability.  The  small  spanwise  irregularity 
in  the  wall  velocity  is  amplified  by  a  factor  7000  which 
is  the  result  of  subcritical  transient  growth  combined 
with  subsequent  growth  due  to  instability.  The  am¬ 
plitude  ratio  at  X  =  500  is  A/Aj  =  654  (figure  15). 
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Figure  17:  Amplitude  of  cross-flow  vortices  at  x  = 
500  in  response  to  forcing  at  different  upstream  posi¬ 
tions,  0  =  —0.4. 

In  general,  an  isolated  or  spanwise  periodic  rough¬ 
ness  element  is  characterized  by  a  Fourier  decomposi¬ 


tion  not  only  in  the  spanwise  direction,  but  also  in  the 
marching  direction  x.  A  more  sophisticated  aiicUysis, 
therefore,  would  allow  for  both  spanwise  and  stream- 
wise  wave  numbers.  Preliminary  results  show  the 
highest  receptivity  for  the  wave  numbers  and  stream- 
wise  position  of  the  critical  disturbance  for  onset  of 
cross-flow  instability. 


3.4  Nonlinear  PSE 

The  previous  section  has  shown  that  the  linear  PSE 
approach  is  an  efficient  and  accurate  method  and  ex¬ 
ceeds  the  capabilities  of  the  local  stability  theory  es¬ 
pecially  with  respect  to  inhomogeneous  problems.  In¬ 
homogeneities  can  not  only  enter  through  boundary 
conditions  but  also  as  a  distributed  force  on  the  right- 
hand  side  of  eqs.  (19).  Such  distributed  forces  most 
naturally  occur  if  we  retain  the  nonlinear  terms  on 
the  right-hand  side  of  eqs.  (7). 

Similar  to  a  spectral  method,  we  assume  the  flow 
field  q'  as  periodic  in  time  t  with  a  period  T  =  2ir/u>i 
and  in  spanwise  direction  z  with  wavelength  = 
27r/^i.  We  expand  the  initial  conditions  and  the  so¬ 
lution  for  X  >  X3  in  a  double  Fourier  series 

q'{x,y,z,t)  - 

00  00 

X:  E  (44) 

71  =  — 00  771=  — CX> 

The  choice  of  ui  and  /?i  is  analogue  to  choosing  the 
computational  domain  in  Navier-Stokes  simulations. 
Some  guidance  is  given  by  controlled  experiments  and 
the  theory  of  primary  and  secondary  instability. 

To  exploit  the  PSE  concept,  we  write  every  mode 
in  terms  of  slowly  varying  functions  of  x, 

CJO 

,  (45a) 

/=0 

,  (456) 

such  that  every  amplitude  function  qnm/(s;,y)  is  as- 
sociated  with  the  wave  function 


Xnml  —  6Xp(2^r;y7<j/ (x)  -b  i^mZ 

where  /?m  —  ctnd  ujn  =  nu)i.  With  the  proper 
indices  attached  to  the  various  quantities,  relations 
(17),  the  linear  left-hand  side  of  eqs.  (18),  and  rela¬ 
tions  (19)  to  (21)  carry  over  to  the  nonlinear  problem 
except  that  eq.  (19)  now  becomes  a  system  of  inho¬ 
mogeneous  equations, 

r  1  ^^.nml 

^nml  qriTn/  i  ^nml 


-h^^y~Nnml  qnm/  =  Ynml  ,  (47) 

ax 

where  the  right-hand  side  Vnmi  =  v) 

sum  over  all  terms  that  contribute  to  the  mode 
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(n,  m,  /).  To  collect  these  terms  from  the  sixfold  sum¬ 
mation  (for  the  quadratic  nonlinearity)  is  straightfor¬ 
ward  with  respect  to  z  and  t,  yet  not  in  x.  The  first 
component  of  Tnmi  is  zero  since  the  continuity  equa¬ 
tion  (7a)  is  linear.  The  terms  that  contribute  to  the 
second  component  of  rnmi  are  obtained  from  eq.  (7b) 
as 


m  —  IJLK,  T 


dut, 


fJ.K  > 


dx 


dUuiiK,  \  .  ^  , 

~TVn  —  u^m  —  ^,X  j  T"  — f/,m  — /x,A  hk\ 


Eqs.  (7c)  and  (7d)  provide  similar  expressions  for  the 
remaining  terms. 

With  respect  to  the  x  direction,  the  decision  is  yet 
open  which  term  contributes  to  which  mode  since 
there  is  no  simple  relation  between  the  indices  A, 
and  /.  To  minimize  the  effect  of  the  PSE  approxima¬ 
tion  and  to  close  the  system  of  equations,  every  mode 
is  subject  to  a  proper  norm,  e.g. 


/„vL,^«.  =  0  (48, 

in  analogy  to  eq.  (34).  This  norm  determines  a^m/ 
and  Onmi  to  make  qnm/  a  slowly  varying  function  oix. 
To  maintain  this  property  in  the  presence  of  Tnmi,  we 
observe  that  the  terms  in  square  brackets  are  slowly 
varying  and  select  /  such  that  for  given  n  and  m 


—  min  — H"  ~  ^nm/|  j  (^9) 

i.e.,  the  term  contributes  to  the  “nearest  neighbor” 
with  respect  to  the  streamwise  wavenumber  to  reduce 
the  residual  variation  with  x.  This  residual  variation 
is  absorbed  into  the  amplitude  function.  Intuitively  it 
is  clear  that  the  effect  of  the  PSE  approximation  can 
be  reduced  by  providing  a  large  number  of  choices  for 
the  selection  (49)  of  /  though  at  higher  computational 
expense. 

In  previous  studies  of  ribbon-controlled  transi¬ 
tion  [19],  successful  runs  were  performed  with  only 
a  single  wavenumber  a^mo  per  mode  and  the  differ¬ 
ence  ~  ^nmo  has  been  neglected. 

However,  this  approach  is  too  limited  for  a  broader 
range  of  applications.  The  need  for  more  wavenum¬ 
bers  anmi  can  be  determined  by  monitoring  A/. 

While  it  is  desirable  to  solve  the  nonlinear  PSE  sys¬ 
tem  in  a  uniform  procedure  for  all  modes,  the  mode 
Qoo  requires  certain  exceptions.  This  special  role  of 
the  mean-flow  distortion  arises  since  it  is  governed 
-  like  the  basic  flow  -  by  the  boundary-layer  equa¬ 
tions,  not  the  Navier-Stokes  equations.  In  general, 
the  wavenumber  associated  with  the  mean-flow  dis¬ 
tortion  is  unnecessary,  ctQoi  =  0.  The  growth  of  the 
distortion  caused  by  all  terms  is  absorbed  into  a  single 
amplitude  function  qooo(^jy)-  The  boundary  condi¬ 
tions  on  the  mean-flow  distortion  are  different  from 


those  of  all  other  modes  since  the  growth  of  distur¬ 
bances  in  the  boundary  layer  causes  changes  in  the 
displacement  thickness  and  hence  t;ooo  ^  0  as  ^  ^  oo. 

In  contrast  to  unsteady  modes  with  n  0,  steady 
modes  (including  the  mean-flow  distortion)  exhibit 
a  pronounced  dependence  on  history  and  accumulate 
small  differences  or  changes  during  the  marching  solu¬ 
tion.  Since  the  mean-flow  and  its  derivatives  strongly 
affect  the  stability  characteristics,  utmost  care  has  to 
be  paid  to  the  treatment  of  the  mean  flow  and  its 
nonlinear  distortion.  Also  note  that  quantities  such 
as  the  skin  friction  coefficient  C/(x)  (and  heat  trans¬ 
fer  coefficient  St{x)  in  compressible  flows)  are  solely 
governed  by  the  mean  flow.  All  periodic  components 
average  to  a  zero  contribution  when  the  local  mean 
is  taken  over  t  and  z. 

3.4.1  Implementation 

For  numerical  applications,  the  Fourier  series  (44)  is 
truncated  at  prescribed  values  of  |n|  <  N,  \m\  <  M 
and  the  physical  flow  field  takes  the  form 

N  M 

q'ix,y,z,t)=  ^  • 

n  =  —  N  Tn—-M 

(50) 

Further,  a  limit  /  <  L  is  specified  to  restrict  the  num¬ 
ber  of  modes  with  different  at  given  n,  m.  While 
the  summation  includes  modes  with  n  <  0,  it  is  un¬ 
necessary  to  compute  these  modes  since  the  physical 
field  must  be  real  and  hence 

=  q'lrn  (51) 

where  f  denotes  the  complex  conjugate.  For  two- 
dimensional  basic  flows  such  as  Blasius  flow,  the  com¬ 
putation  can  be  further  reduced  by  the  assumption  of 
symmetry  in  z, 

^  n,  — m  —  [^nrTD^nrri}  ^nmjPnm] 

This  symmetry  can  be  derived  from  the  governing 
equations  but  may  be  used  only  if  initial  and  bound¬ 
ary  conditions  are  symmetric.  Without  assuming 
symmetry,  the  PSE  code  can  track  the  evolution  of 
asymmetric  initial  conditions,  e.g.  of  a  single  oblique 
wave.  For  three-dimensional  boundary  layers,  eq. 
(52)  is  invalid  if  both  oj  ^  0  and  0^0. 

The  “best  approach”  to  the  treatment  of  the  non¬ 
linear  terms  during  the  marching  procedure  requires 
the  usual  compromise  between  accuracy  and  effi¬ 
ciency.  For  many  calculations  it  is  sufficient  to  treat 
the  terms  explicitly.  The  explicit  treatment  is  more 
efficient  since  the  computation  of  the  nonlinear  in¬ 
teraction  terms  can  be  time  consuming.  To  increase 
the  accuracy  of  the  results,  the  marching  step  for  the 
explicit  integration  should  be  small.  Small  march¬ 
ing  steps  may  require  suppression  of  the  term  pj,  to 
remove  the  small  upstream  influence.  For  more  ac¬ 
curate  calculations  and  studies  far  into  the  transition 
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Figure  18:  Schematic  of  the  modes.  The  modes  in 
the  dashed  frame  constitute  the  physical  field.  The 
dark  shaded  area  is  excluded  from  the  computation. 
The  modes  in  the  light  shaded  area  can  be  generated 
by  symmetry  in  two-dimensional  flows. 

regime,  the  implicit  treatment  is  preferable  although 
the  cost  per  iteration  step  is  higher.  Often  it  is  neces¬ 
sary  to  recalculate  the  nonlinear  terms  and  the  mean- 
flow  distortion  after  updating  the  wavenumber  for  ev¬ 
ery  single  mode  to  achieve  convergence  for  strongly 
growing  modes.  The  increased  effort  per  iteration 
step  is  often  outweighed  by  more  rapid  convergence 
and  the  opportunity  to  use  larger  marching  steps. 

Initial  data  are  usually  given  only  for  a  few  modes, 
e.g.  for  modes  (2,0)  and  (1,0)  as  indicated  by  the 
dark  circles  in  figure  18.  The  other  modes  (shaded 
circles)  are  generated  as  the  nonlinear  forcing  terms 
come  into  play.  To  avoid  strong  transients  for  initial 
data  with  large  amplitudes,  it  is  useful  to  calculate 
the  nonlinear  terms  and  generate  at  every  position 
the  yet  missing  modes  using  the  inhomogeneous  ver¬ 
sion  of  the  local  stability  equations  if  the  forcing  ex¬ 
ceeds  a  low  threshold.  This  procedure  provides  rea¬ 
sonable  starting  values  for  the  first  wavenumber  up¬ 
date.  For  the  mean-flow  distortion,  this  approach  is 
inapplicable.  Various  methods  of  generating  the  dis¬ 
torted  mean  flow  at  the  starting  position  have  been 
tested  and  compared.  Since  the  mean  flow  (and  other 
steady  modes)  exhibit  a  pronounced  memory,  differ¬ 
ent  choices  at  the  initial  position  have  a  lasting  effect 
on  the  nonlinear  evolution  of  the  flow  which  may  be 
significant  in  runs  with  large  initial  amplitudes.  It 
appears  best  to  assume  the  mean  flow  distortion  to 
be  zero  upstream  of  the  initial  position.  This  choice 
is  well  justified  for  runs  with  small  initial  amplitudes. 

Many  applications  such  as  the  example  in  figure  18 
leave  certain  modes  unused  (white  circles).  In  the  in¬ 
terest  of  studying  advanced  problems  with  limited  re¬ 
sources,  this  fact  can  be  exploited  by  a  flexible  mem¬ 
ory  allocation.  For  research  applications,  it  is  use¬ 
ful  to  associate  the  modes  with  two  tables,  the  first 
of  which  keeps  track  of  existence  and  status  of  each 
mode  while  the  second  table  indicates  pre-determined 
properties.  By  assigning  special  properties,  modes 
can  receive  special  treatment  or  be  prohibited  to  par¬ 


ticipate  in  nonlinear  interactions. 

Most  results  shown  in  the  following  have  been  pro¬ 
duced  as  tutorial  examples  to  illustrate  capabilities 
and  phenomena.  No  special  care  has  been  applied  to 
produce  archival  results.  Initial  conditions  have  been 
generated  by  the  modified  LST  with  the  transverse 
velocity  terms  included. 

3.4.2  Nonlinear  Effects  on  Single  TS  Modes 

The  nonlinear  evolution  of  single  modes  was  once  the 
focus  of  weakly  nonlinear  theories  which  were  limited 
by  the  inability  to  incorporate  both  nonlinearity  and 
nonparallelism  [53,  55].  This  inability  especially  af¬ 
fected  the  mean-flow  distortion  and  prevented  higher- 
order  perturbation  analysis.  The  nonlinear  evolution 
of  a  TS  wave  wcls  also  the  corner  stone  of  spatial 
Navier-Stokes  solutions  [56].  Therefore,  this  case  has 
also  been  chosen  as  a  first  application  of  the  PSE 
approach  [19,  57,  22].  The  PSE  results  for  a  TS 
wave  with  lji  =  0.0344  and  an  initial  rms  ampli¬ 
tude  of  0.25%  at  Reo  =  400  were  found  in  excellent 
agreement  with  DNS  results.  Results  for  the  same 
case  were  also  presented  by  Chang  et  al.  [21]  and 
also  found  in  agreement  with  DNS  results  although 
their  maximum  amplitude  is  somewhat  higher  than 
the  value  reported  by  Bertolotti  et  al.  [22]. 
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Figure  19:  Linear  and  nonlinear  evolution  of  a  TS 
wave  of  frequency  wi  =  0.0344  with  =  0.0025 
at  Xs  —  Reo  =  400.  The  various  modes  are  labeled 
n,  m  for  their  frequency  u>n  and  spanwise  wavenumber 
Prn  as  in  figure  18  (/  is  omitted  if  L  =  0). 

Figure  19  shows  a  similar  result  obtained  with 
iV  =  5  for  the  same  case  in  comparison  with  the  lin¬ 
ear  result.  Although  nonlinearity  is  small,  as  shown 
by  the  low  amplitudes  of  the  harmonics,  it  causes  a 
significant  change  in  the  peak  amplitude  and  the  lo¬ 
cation  of  branch  II.  The  destabilization  near  branch 
II  is  consistent  with  the  results  of  the  weakly  nonlin¬ 
ear  theory  [53,  54],  direct  numerical  simulation  [58], 
and  asymptotic  theory  [59].  Notable  is  the  evolution 
of  the  mean-flow  distortion,  (0,0).  The  “kink”  near 
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X  —  1700  is  an  artifact  of  plots  of  the  maximum  abso¬ 
lute  value  for  oscillatory  functions.  Different  parts  of 
the  curve  are  associated  with  peaks  at  different  dis¬ 
tance  from  the  plate.  While  plots  of  local  amplitudes 
are  often  necessary  to  compare  with  measurements, 
plots  based  on  the  amplitude  A  or  the  kinetic  energy 
are  in  general  easier  to  interpret. 

The  mean-flow  distortion  assumes  a  maximum  fur¬ 
ther  downstream  than  the  nonlinear  TS  wave.  In 
other  words,  the  mean-flow  distortion  continues  to 
accumulate  the  nonlinear  effects  in  the  presence  of  a 
finite-amplitude  wave  even  if  this  wave  decays.  While 
traveling  waves  quickly  adapt  to  local  conditions,  the 
mean-flow  distortion  and  other  steady  modes  possess 
a  pronounced  memory.  This  memory  cannot  be  mod¬ 
elled  by  the  local  weakly  nonlinear  analysis  but  re¬ 
quires  solving  an  initial-value  problem. 

The  maximum  amplitude  of  the  TS  wave  in  fig¬ 
ure  19  is  slightly  lower  than  the  result  of  Bertolotti 
et  al.  [22].  This  difference  is  largely  due  to  the  use 
of  different  initial  conditions;  we  use  a  linear  stability 
mode  while  results  of  the  local  weakly  nonlinear  or 
nonparallel  theories  were  used  in  earlier  work.  The 
sensitivity  of  the  solution  to  initial  conditions  is  a 
fundamental  property  of  the  stability  and  transition 
prob  ems  we  deal  with.  Every  small  change  in  ini¬ 
tial  conditions  or  numerical  treatment  causes  minute 
changes  in  the  growth  rates  which  accumulate  over 
the  long  distance  downstream  to  small  yet  notice¬ 
able  differences  in  the  amplitude.  Different  initial 
conditions  also  have  caused  the  different  results  of 
Bertolotti  [19]  and  Chang  et  al.  [21]  (M.  R.  Malik, 
personal  communication,  1993). 

A  second  principal  difference  between  the  earlier 
and  present  approach  is  the  update  of  the  wave  num¬ 
ber.  Leaning  on  the  weakly  nonlinear  theory  for  par¬ 
allel  flow  [60],  the  earlier  formulation  implemented 
the  restriction 

(Q^r)nO  ^  ^(^r)lO 

on  the  wavenumber  and  updated  only  the  growth  rate 
from  an  integral  norm.  In  nonparallel  flows  there 
is  no  rigorous  justification  for  this  property.  There¬ 
fore,  our  analysis  determines  both  wavenumber  and 
growth  rate  by  the  norm.  In  the  case  shown  above, 
relation  (53)  is  satisfied  only  within  about  1%  once 
the  initial  transients  decayed.  Although  this  devia¬ 
tion  may  be  negligible  and  the  use  of  eq.  (53)  saves 
the  cost  for  iterating  on  the  higher  wavenumbers,  this 
approach  is  not  generally  applicable  and  the  assump¬ 
tion  of  a  single  dominating  (“ribbon-induced”)  mode 
is  too  restrictive. 

3.4.3  Secondary  Instabilities 

The  next  major  test  for  the  PSE  approach  and  its  im¬ 
plementation  is  the  ability  to  track  the  simultaneous 
evolution  of  multiple  modes  and  their  interactions, 


especially  for  the  parametric  resonance  mechanisms 
that  lead  to  secondary  instabilities  [4]. 

Earlier  PSE  studies  [19,  57]  considered  the  effect  of 
a  finite-amplitude  TS  wave  on  a  subharmonic  mode 
of  secondary  instability,  i.e.  the  initial  conditions  for 
the  subharmonic  mode  were  obtained  from  the  Flo- 
quet  analysis  of  secondary  instability.  Suppressing 
the  feedback  of  the  subharmonic  mode  on  the  TS 
wave,  the  subharmonic  mode  evolved  in  close  prox¬ 
imity  to  the  predictions  of  the  Floquet  theory. 


Figure  20:  Linear  subharmonic  secondary  instability 
and  influence  of  the  initial  phase. 

Figure  20  shows  a  similar  result  for  the  conditions 
of  the  experiment  of  Kachanov  h  Levchenko  [10], 
X,  =  /?eo  =  420,  =  0.005208,  A  =  0.14.  Here, 

however,  the  initial  data  are  given  by  the  two  pri¬ 
mary  modes  (10,  0)  and  (5, 1)  with  rms  amplitudes  of 
0.46%  and  0.0035%,  respectively.  Two  cases  with  dif¬ 
ferent  phase  shifts  between  the  two  modes  are  shown. 
The  feedback  of  the  three-dimensional  wave  on  the 
two-dimensional  mode  and  the  mean-flow  distortion 
are  suppressed. 

In  absence  of  the  TS  wave,  the  primary  mode  (5,1) 
would  decay  as  shown.  Affected  by  the  TS  wave,  the 
three-dimensional  mode  exhibits  a  strong  transient, 
and  the  wavenumber  settles  near  twice  the  wavenum¬ 
ber  of  the  TS  wave  while  the  growth  rate  strongly 
increases.  The  perfect  synchronization  of  both  waves 
obtained  by  Floquet  analysis  (and  enforced  in  ear¬ 
lier  work)  is  not  achieved  in  the  nonparallel  flow 
with  streamwise  varying  TS  amplitude.  The  near- 
synchronization  is  lost  at  X  ^  2000  and  the  subhar¬ 
monic  mode  develops  on  its  own  while  the  parametric 
excitation  has  faded  away.  The  subharmonic  growth 
rates  are  slightly  smaller  than  the  theoretical  predic¬ 
tions.  The  transient  evolution  and  the  maximum  am¬ 
plitude  of  the  subharmonic  mode  depend  on  the  phase 
relation  to  the  TS  wave. 

Closely  related  to  subharmonic  secondary  instabil¬ 
ity  is  the  mechanism  of  combination  resonance  in  the 
presence  of  the  TS  wave  (10,0)  and  a  detuned  mode 
(4,  1)  with  amplitudes  and  other  parameters  as  above, 
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and  the  same  phase.  As  figure  21  shows,  the  one¬ 
sided  nonlinear  interaction  introduces  the  mode  (6, 1) 
such  that  the  sum  of  the  frequencies  of  the  three- 
dimensional  modes  equals  the  frequency  of  the  TS 
wave.  For  700  <  x  <  1300  the  combination  modes  de¬ 
velop  with  similar  large  growth  rates.  Further  down¬ 
stream,  where  the  rms  amplitude  of  the  TS  wave  is 
less  than  0.3%,  each  of  the  three-dimensional  modes 
is  governed  by  its  own  linear  stability  characteristics. 
The  maximum  amplitudes  due  to  combination  res¬ 
onance  is  less  than  for  subharmonic  resonance.  A 
similar  result  is  obtained  by  specifying  mode  (6, 1) 
instead  of  (4, 1)  by  the  initial  conditions. 


Figure  21:  Combination  resonance. 

The  third  type  of  secondary  instability,  fundamen¬ 
tal  or  principal  parametric  resonance  is  rather  weak 
at  the  relatively  low  TS  amplitudes  in  the  Kachanov- 
Levchenko  experiment.  We,  therefore,  consider  this 
type  at  lower  frequency  where  the  TS  waves  experi¬ 
ence  stronger  growth.  Most  of  the  microscopic  exper¬ 
iments  on  K-type  breakdown  [8,  64]  were  performed 
at  these  lower  frequencies. 

The  K-type  of  secondary  instability  is  caused  by 
the  interaction  of  a  TS  wave  (1,0)  with  either  a  pair 
of  oblique  waves  (1,1)  or  a  longitudinal  vortex  (0,1). 
The  latter  interaction  is  favored  by  KlebanofF  modes 
(caused  by  “spacers”  underneath  the  ribbon  in  exper¬ 
iments). 

Figure  22  shows  the  evolution  of  the  amplitudes 
from  initial  modes  (1,0),  (1,1)  with  Xs  —  Reo  =  600, 
0^1  =  0.03,  /3i  =  0.15,  and  small  initial  rms  ampli¬ 
tudes  of  0.1%  and  0.01%,  respectively.  The  interac¬ 
tion  generates  mode  (0,1),  a  steady,  spanwise  periodic 
mode.  Because  of  the  different  streamwise  wavenum¬ 
bers  of  the  initial  modes,  the  amplitude  oscillates  and 
increases  in  the  mean  as  the  TS  amplitude  grows.  As 
modes  (0,1)  and  (1,1)  reach  comparable  magnitude, 
they  engage  in  an  interaction  with  each  other  and 
with  the  TS  wave.  Mode  (1,1)  synchronizes  with  the 
TS  wave  and  the  oscillations  of  mode  (0,1)  disappear, 
leaving  a  longitudinal  vortex  mode.  The  spanwise  pe¬ 
riodic  disturbance  grows  rapidly  (a,-  «  0.035)  to  large 


Figure  22:  Fundamental  (K-type)  secondary  instabil¬ 
ity  originating  from  initial  modes  (1,0)  and  (1,1). 

amplitude.  The  unrealistic  amplitudes  in  these  runs 
appear  because  the  effects  of  the  secondary  distur¬ 
bances  on  the  TS  wave  and  the  generation  of  har¬ 
monics  are  suppressed. 

3.4.4  Nonlinear  Receptivity  Studies 

Another  application  of  the  nonlinear  PSE  is  the 
study  of  certain  receptivity  mechanisms  that  in¬ 
volves  the  interaction  of  small  disturbances  at  the 
boundaries,  such  as  the  interaction  of  sound  with  a 
wavy  wall.  Based  on  the  local  theory,  these  mecha¬ 
nisms  have  been  analyzed  by  Crouch  [61]  and  Choud- 
hari  &  Street  [62]  for  two-dimensional  boundary  lay¬ 
ers.  Nonlinearity  is  necessary  to  combine  the  Stokes 
wave  of  frequency  uj  and  the  steady  component  of 
wavenumber  d  induced  by  the  waviness  into  a  forc¬ 
ing  term  for  TS  waves. 

Local  and  PSE  results  have  been  compared  by 
Crouch  &  Bertolotti  [63]  and  found  in  good  agree¬ 
ment.  Small  differences  are  attributed  to  non¬ 
parallelism.  In  these  PSE  computations,  the  forced 
TS  wave  maintains  the  wavenumber  d  and  the  grow¬ 
ing  difference  between  d  and  the  wavenumber  ars 
of  a  TS  eigenmode  is  absorbed  into  the  amplitude 
function. 

Analysis  of  the  simple  forced  linear  oscillator  dis¬ 
cussed  by  Crouch  [61]  shows,  that  d  dominates  the 
behaviour  of  the  forced  mode  only  in  the  stable  do¬ 
main.  After  passing  through  resonance  (at  or  near  the 
neutral  point),  the  wavenumber  of  the  forced  mode 
approaches  that  of  the  eigenmode.  Therefore,  we  pre¬ 
fer  to  update  the  wavenumber  of  the  forced  compo¬ 
nent  and  extract  the  oscillatory  behaviour  from  the 
amplitude  function. 

The  results  in  figure  23  are  for  Rcq  =  400,  x,  = 
300,  u  =  0.0224  (F  =  56  •  10”®)  and  a  wall  waviness 
with  d  =  0.0696944  (as  in  [63]).  The  difference  be¬ 
tween  the  two  evolutions  with  and  without  updating 
the  wavenumber  of  the  forced  TS  wave  is  consider- 
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Figure  23:  Interaction  of  sound  with  a  wavy  wall. 
Effect  of  the  wavenumber  update. 

able  and  cannot  be  attributed  to  different  effects  of 
the  PSE  approximation.  At  small  Reynolds  numbers 
up  to  Re  700,  the  smaller  response  with  wavenum¬ 
ber  update  is  due  to  the  changing  phase  difference  be¬ 
tween  the  forced  wave  and  the  forcing  term  which  also 
causes  an  ondulation  of  the  amplitude  that  diminishes 
cLs  the  natural  oscillation  takes  over  and  a  — C(ts- 
Although  the  difference  in  the  wavenumbers  a  and 
ar  is  small,  it  causes  different  growth  characteristics 
at  higher  Re  as  shown  by  the  shift  in  the  branch  II 
position. 

3.5  Ribbon-Controlled  Transition 

A  more  severe  test  of  the  nonlinear  PSE  is  the  repro¬ 
duction  of  the  microscopic  experiments  of  Klebanoff, 
Tidstrom  k  Sargent  [8],  Kachanov  k  Lechenko  [10], 
Cornelius  [64],  and  others  on  transition  in  the  flat- 
plate  boundary  layer.  These  experiments  are  also  the 
benchmark  cases  for  spatial  DNS  codes  and  serve  to 
analyze  the  nonlinear  development  of  transition  from 
different  types  of  secondary  instability.  Besides  for 
these  test  cases,  the  PSE  can  be  used  to  study  re¬ 
lated  problems,  in  particular  problems  that  involve 
numerous  modes  such  as  wave  packets. 

PSE  studies  for  these  test  cases  have  been  per¬ 
formed  by  Bertolotti  [19,  57]  and  compared  with  ex¬ 
periments  and  Navier-Stokes  solutions.  Good  agree¬ 
ment  has  been  found  except  for  some  experimental 
data  as  the  2/  mode  (4,0)  of  Kachanov  k  Levchenko. 
For  this  mode,  both  PSE  and  DNS  show  a  continuous 
rise  while  the  experiment  found  a  significant  decline 
followed  by  a  steep  rise  near  onset  of  transition.  The 
discrepancy  was  attributed  to  experimental  error. 

Of  the  various  analyses  of  transition  arising  from 
H-type,  K-type,  or  combination  resonance  or  interac¬ 
tion  of  two  oblique  modes  (1,1)  and  (1,-1),  we  show 
only  one  example^  that  helps  to  clarify  the  earlier 
discrepancy  between  the  experiment  of  Kachanov  k 
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Levchenko  [10]  and  computations  [65,  19,  57], 

PSE  studies  of  this  case  by  Bertolotti  used  a  dom¬ 
inant  TS  wave  and  a  subharmonic  secondary  mode 
from  Floquet  analysis  as  initial  conditions.  Both 
modes  were  in  the  correct  phase  and  synchronization. 
Here  we  continue  our  study  of  subharmonic  secondary 
instability  (section  3.4.3)  by  fully  accounting  for  non¬ 
linearity.  Deviating  from  previous  work,  we  use  a  pri¬ 
mary  mode  to  introduce  a  subharmonic  component. 
This  choice  allows  tracing  the  synchronization  pro¬ 
cess  between  two-dimensional  and  three-dimensional 
components  and  the  effect  of  the  phase  difference  be¬ 
tween  the  initial  modes.  All  wavenumbers  are  de¬ 
termined  by  the  norm  without  any  assumption  on 
synchronization. 


Figure  24:  Amplitude  growth  of  selected  modes  in  a 
PSE  simulation  of  the  Kachanov-Levchenko  experi¬ 
ment. 

Figure  24  shows  the  nonlinear  evolution  {N  = 
4,  M  =  2)  of  selected  modes  from  the  same  initial 
conditions  as  for  figure  20  (phase  0'’).  With  full  in¬ 
teraction  permitted,  the  growing  subharmonic  mode 
(1,1)  prevents  the  decline  of  the  TS  wave.  Here  and 
in  other  runs  with  different  truncations  and  initial 
phase,  mode  (4,0)  behaves  as  in  the  experiment.  The 
earlier  disagreement  between  PSE  results  and  experi¬ 
mental  data  is  apparently  due  to  the  particular  choice 
of  initial  amplitudes  and  the  enforced  resonance  from 
the  outset. 

3.6  Nonlinear  Vortices 

Transition  evolving  from  Gortler  vortices  or  cross- 
flow  vortices  is  still  a  matter  of  experimental  and 
computational  research.  These  types  of  disturbances 
do  not  initiate  secondary  instabilities  at  small  ampli¬ 
tudes  where  the  approximations  of  the  Floquet  anal¬ 
ysis  [4]  are  justified.  Instead,  vortices  grow  to  large 
amplitudes  and  develop  strong  harmonics  before  sec¬ 
ondary  instabilities  develop  on  the  embedded  layers 
of  high  shear. 


4-23 


I 

I 


Figure  25:  Nonlinear  growth  of  a  steady  cross-flow 
vortex  in  swept  Hiemenz  flow. 

The  nonlinear  evolution  of  a  steady  cross-flow  vor¬ 
tex  in  swept  Hiemenz  flow  for  the  conditions  of  fig¬ 
ure  15,  Re  =  500,  /?  =  —0.4,  x,  =  190,  an  initial 
amplitude  of  \u\max  =  10~^/>/2,  and  M  =  8  is  shown 
in  figure  25.  The  rapid  growth  of  mode  (0,1)  causes 
the  rise  of  harmonics  until  saturation  is  reached  near 
a?  600.  For  x  >  700,  the  maximum  of  the  mean-flow 
distortion  |uoo|  exceeds  10%  and  the  spanwise  peri¬ 
odic  component  |iioi|  exceeds  17%  of  the  reference 
velocity  We  ♦  These  ratios  are  smaller  with  respect  to 
Ue  =  xWe! Re.  saturation  level  is  insensitive  to 

the  initial  amplitude.  For  higher  initial  amplitudes, 
saturation  is  reached  at  smaller  x.  Similar  character¬ 
istics  have  been  found  for  unsteady  cross-flow  vortices 
(see  also  [46,  47]). 


Figure  26:  Contour  lines  of  equal  u  velocity  at  ar  = 
716.  The  contour  levels  are  multiples  of  5%  of  the 
edge  velocity  Ue- 

Figure  26  shows  a  cross  section  of  the  vortex  at 
X  =  716  over  one  period  in  z  and  y  <  8  (to  scale). 
The  different  spacing  of  the  contour  lines  near  the 
wall  indicates  a  strong  spanwise  variation  of  the  wall- 
shear  stress.  A  vertical  cut  near  the  center  reveals  a 
strongly  inflectional  velocity  distribution. 

While  the  general  features  of  figure  25  can  be  ob¬ 
tained  with  M  =  2,  the  detailed  structure  of  the 
cross-flow  vortex  after  saturation  depends  on  the 


number  of  harmonics  taken  into  account.  The  rel¬ 
atively  strong  influence  of  harmonics  raises  questions 
on  the  validity  of  secondary  instability  results  ob¬ 
tained  at  low  truncation.  Floquet  analysis  of  the  flow 
at  fixed  X  under  neglect  of  the  streamwise  variation 
but  accounting  for  all  relevant  harmonics  is  a  com¬ 
putationally  expensive  task.  Alternatively,  the  lin¬ 
ear  primary  instability  of  the  three-dimensional  flow 
at  fixed  X  can  be  analyzed  using  two-dimensional 
eigenfunctions  in  z  and  y,  but  the  computational  de¬ 
mand  remains  high.  At  this  time,  reliable  results  on 
the  instability  of  highly  nonlinear  cross-flow  vortices 
are  unavailable.  As  a  consequence,  the  breakdown 
of  these  vortices  and  the  resulting  transition  process 
have  not  yet  been  analyzed.  Work  toward  closing  this 
gap  is  in  progress. 

4  PSE  for  Compressible  Flows 

The  successful  application  of  the  PSE  for  incompress¬ 
ible  flows  and  the  demand  for  improved  predictive 
tools  in  aerodynamic  design  suggest  the  extension 
to  compressible  flows.  While  this  extension  appears 
straightforward,  it  faces  a  series  of  hurdles  and  deci¬ 
sions: 

•  The  sheer  size  of  the  nonlinear  stability  equa¬ 
tions  and  associated  code  that  need  careful 
checks  and  support  by  symbolic  manipulation 
even  in  the  simplest  case  of  Cartesian  coordi¬ 
nates. 

•  The  implementation  of  general  curvilinear  coor¬ 
dinates  to  adapt  to  realistic  geometries. 

•  The  choice  of  the  “best”  set  of  variables  between 
p,  T,  or  p,  primitive  or  conservative  variables, 
covariant  or  contravariant  velocity  components. 

•  The  choice  of  models  for  the  thermodynamic 
properties  of  the  working  gas  at  high  Mach  num¬ 
bers. 

•  The  failure  or  inefficiency  of  proven  numerical 
methods. 

•  The  uncertain  accuracy  of  the  basic  flow  and 
its  derivatives  generated  by  traditional  compu¬ 
tational  methods  and  bound  ary- layer  codes  for 
flows  over  aerodynamic  bodies. 

•  The  need  for  consistent  interpolation  methods  to 
obtain  the  basic-flow  data  on  the  grid  used  for 
stability  analysis  from  the  data  given  on  some 
other  computational  grid. 

•  The  unavailability  of  detailed  test  data  from  mi¬ 
croscopic  experiments  or  spatial  DNS  at  high 
speeds. 

•  The  lack  of  physical  understanding  for  mecha¬ 
nisms  and  results. 
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As  a  rule  of  thumb,  the  difficulties  increase  with 
the  Mach  number  and  the  complexity  of  geometry 
and  disturbance  environment. 

4.1  Stability  Equations 

The  PSE  for  compressible  flow  in  Cartesian  coordi¬ 
nates  up  to  quadratic  nonlinear  terms  were  devel¬ 
oped  by  Bertolotti  and  applied  to  study  the  non¬ 
parallel  linear  stability  of  the  flat-plate  boundary 
layer  [19,  66,  67].  Based  on  his  work,  the  PSE  ap¬ 
proach  was  also  adopted  by  Chang  et  al.  [21]  to  per¬ 
form  a  nonlinear  analysis  of  this  flow.  While  the 
choice  of  variables  and  coordinates  for  the  traditional 
stability  analysis  is  important  because  of  the  effect 
on  the  parallel-flow  assumption,  it  is  less  critical  for 
the  linear  PSE.  However,  the  use  of  p  and  T  is  ad¬ 
vantageous  for  the  nonlinear  formulation.  The  com¬ 
plete  nonlinear  formulation  including  cubic  and  quar- 
tic  nonlinearities  in  general  curvilinear  coordinates 
has  been  developed  by  Stuckert  et  al.  [68].  The  ve¬ 
locity  field  is  described  by  the  contravariant  velocity 
components. 

Our  implementation  of  the  nonlinear  compressible 
PSE  is  rather  general.  Different  sets  of  disturbance 
variables,  e.g.  q  =  [p,  T,  (pu)^]^  or  q  =  [p,T,  u^]^ 
can  be  chosen  for  the  analysis.  Different  models  of 
the  thermodynamic  properties  can  be  chosen,  which 
is  helpful  in  comparing  with  previous  work.  Two 
separate  time-dependent  coordinate  transformations 
y,  z,  f)  ,  k  =  1,2,  3  can  be  specified:  for  the 
computational  grid  of  the  basic  flow  and  for  the  grid 
used  in  the  PSE  analysis.  This  capability  enables  use 
of  the  same  core  modules  for  diverse  problems  such  as 
flat  plate  or  rotating  disk,  swept  wings  or  blunt  cones. 
Grouped  around  this  core  are  application-dependent 
modules  as  well  as  interfaces  to  similarity  solutions, 
boundary-layer  codes,  and  the  PlotSD  files  produced 
by  most  computations  of  the  basic  flow. 


4.2  Numerical  Aspects 

While  the  incompressible  problem  spends  most  of  the 
computer  time  in  solving  linear  algebraic  systems, 
the  setup  of  these  systems  for  nonlinear  compressible 
analysis  consumes  a  major  fraction  of  the  total  time. 
Obviously,  the  setup  time  increases  with  the  number 
of  grid  points  needed  to  resolve  the  flow  in  normal 
to  the  boundary  and  with  the  order  of  approximation 
(quadratic,  cubic,  or  quartic)  to  the  nonlinear  terms. 

The  reliable  and  accurate  single-domain  spectral 
method  successfully  used  for  the  incompressible  for¬ 
mulation  in  terms  of  Orr-Sommerfeld  and  Squire 
equation  [19]  loses  advantage  for  the  compressible 
problem  because  the  dimension  of  the  algebraic  sys¬ 
tem  with  a  completely  filled  matrix  increases  from 
2J  to  57  where  J  is  the  number  of  collocation  (or 
grid)  points  normal  to  the  boundary.  At  high  Mach 


numbers,  the  inflexible  arrangement  of  the  colloca¬ 
tion  points  near  the  boundary  makes  it  difficult  to 
resolve  details  of  higher  modes  with  critical  layers 
near  the  edge  of  the  thick  boundary  layers.  Bertolotti 
(personal  communication,  1993)  has  adopted  a  multi- 
domain  spectral  method  to  overcome  these  shortcom¬ 
ings.  While  the  multi-domain  approach  reduces  the 
matrix  size,  it  demands  additional  numerical  speci¬ 
fications  which  need  optimization  and  distract  from 
the  physical  problem. 

The  fourth-order  compact  (Euler-Maclaurin) 
method  favored  by  Chang  et  al.  [21]  requires  the  te¬ 
dious  reformulation  of  the  PSE  equations  as  a  system 
of  first-order  differential  equations  in  f  ^  and  leads  to 
numerous  terms  involving  the  inverse  (I/^)~^  of  the 
small  transverse  component  of  the  basic  flow.  In 
computational  basic  flows,  this  component  is  usually 
inaccurate.  Moreover,  can  change  sign  and  lead 
to  unpleasant  divisions  by  zero. 

We  maintain  the  single-domain  spectral  method 
for  special-purpose  studies  such  as  obtaining  spec¬ 
tra  of  eigenvalues  or  highly  accurate  results  for  com¬ 
parisons.  The  routine  calculations  are  based  on 
Hermitian  finite-difference  methods  [70]  for  second- 
order  differential  equations.  For  nonuniform  grids, 
two-point  formulas  permit  overall  accuracy  up  to 
fourth  order,  three-point  formulas  up  to  seventh  or¬ 
der.  Nonuniform  grids  with  moderate  stretching  re¬ 
duce  the  number  of  grid  points  needed  to  achieve  a 
certain  accuracy  at  lower  Mach  numbers.  Equidistant 
grids  are  preferable  at  high  Mach  numbers, 
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Figure  27:  Error  of  the  temporal  eigenvalue  calcu¬ 
lated  with  the  three-point  Hermitian  method  at  high 
Mach  number.  The  dashed  line  indicates  seventh- 
order  accuracy. 

Figure  27  shows  the  result  of  a  convergence  study 
for  test  case  IV  of  Malik  [69]  at  Ma  =  10,  Re  =  1000, 
Te  =  lll.lll^A",  and  a  =  0.12  for  the  compress¬ 
ible  flat-plate  boundary  layer.  Our  “exact”  result  for 
the  temporal  eigenvalue  was  obtained  by  the  spec¬ 
tral  method  with  J  >  50  as  w  =  0.115861436  + 


20.000137657  (different  from  the  spectral  results  given 
in  ref.  [69]).  The  error  comprises  both  the  errors  of 
the  basic-flow  calculation  and  stability  analysis.  To 
obtain  the  correct  eigenvalue  within  an  absolute  er¬ 
ror  of  10“^  or  10“®  requires  about  100  or  200  grid 
points,  respectively.  The  spectral  result  is  matched 
with  300  points.  At  lower  Mach  numbers,  Ma  <  2, 
good  accuracy  can  be  achieved  with  as  few  as  50  grid 
points. 

We  have  compared  application  of  the  Hermitian 
formulas  with  and  without  elimination,  and  with 
block-tridiagonal  or  banded  matrix  solvers  without 
finding  a  superior  approach.  Using  banded  matrix 
solvers  without  elimination  requires  more  time  to 
solve  the  system  but  saves  the  time  for  the  elimi¬ 
nation  step  during  setup  and  for  retrieving  the  sec¬ 
ond  derivatives  needed  when  calculating  the  nonlin¬ 
ear  terms.  A  clear  advantage,  however,  is  the  use  of 
the  LAPACK  library  subroutines  and  of  the  BIAS 
libraries  which  are  optimized  for  various  supercom¬ 
puter  and  workstation  architectures. 

While  we  have  overcome  the  numerical  difficulties 
in  solving  the  stability  problem  for  exact  or  similar¬ 
ity  solutions  as  basic  flows,  applications  to  computa¬ 
tional  solutions  for  realistic  geometries  face  additional 
numerical  problems.  These  problems  often  start  with 
insufficient  accuracy  of  the  coordinates  that  define 
the  shape  of  the  body,  e.g.  an  airfoil.  Small  oscil¬ 
lations  in  the  pressure  gradient  and  the  curvature 
along  the  body  cause  strong  oscillations  in  the  sta¬ 
bility  characteristics. 

Most  computational  methods  for  aerodynamic  ap¬ 
plications  are  second-order  accurate  and  designed  to 
provide  reliable  results  for  the  pressure  distribution 
at  the  surface.  The  stability  characteristics,  however, 
are  governed  by  the  whole  flow  field  in  the  boundary 
layer  and  beyond.  Our  experience  shows  that  sat¬ 
isfactory  accuracy  can  be  obtained  with  boundary- 
layer  codes.  Viscous  flow  solutions  from  Navier- 
Stokes  or  PNS  solvers  often  lack  the  accuracy  needed 
for  a  reliable  stability  analysis. 

4.3  Validation  of  the  PSE  Code 

To  guarantee  proper  function  of  the  compressible 
PSE  code,  extensive  tests  have  been  performed  for 
the  speed  range  from  Ma  0  to  Ma  —  10.  These 
tests  range  from  comparison  with  independent  results 
to  solving  the  same  problem  in  different  coordinate 
systems  and  with  different  sets  of  variables. 

Figure  28  shows  the  comparison  with  results  of  El- 
Hady  [71]  for  subharmonic  secondary  instability  in 
the  flat-plate  boundary  layer  at  Ma  =1.2.  Two  PSE 
runs  were  performed  using  Cartesian  and  similarity 
coordinates,  respectively.  Although  the  PSE  runs 
have  been  tailored  to  match  the  locally-parallel  flow 
assumption,  the  results  deviate  from  El-Hady^s  due 
to  a  slight  shift  of  the  instability  region  for  the  pri¬ 
mary  mode  toward  higher  Reynolds  number  R.  Ex- 
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Figure  28:  Subharmonic  secondary  instability  at 
Ma  =  1.2.  Comparison  of  the  PSE  solution  in  Carte¬ 
sian  (solid)  and  similarity  coordinates  (dashed)  with 
results  ofEl-Hady  (symbols).  □:  PSE,  primary  mode 
(2,0),  C>:  PSE,  subharmonic  mode  (1,1).  x:  El-Hady, 
primary  mode  (2,0),  *:  El-Hady,  subharmonic  mode 
(1.1). 

cept  for  the  initial  transients,  the  PSE  results  agree. 
The  transients  are  stronger  in  Cartesian  coordinates 
because  the  initial  data  are  more  affected  by  the 
parallel-flow  assumption  of  the  LST. 

5  Practical  Applications 

While  most  of  our  work  on  compressible  flows  has 
served  to  develop  and  implement  the  PSE  approach 
together  with  efficient  and  accurate  numerical  meth¬ 
ods  for  integration  and  interpolation,  to  validate  the 
code,  and  gain  an  overview  of  the  differences  between 
LST  and  PSE  results,  we  have  also  performed  feasi¬ 
bility  studies  and  tested  the  code  in  some  real-world 
applications. 

5.1  Flow  over  Swept  Wings 

Flows  over  the  swept  wings  of  airplanes  are  routinely 
analyzed  for  the  transition  location  using  the  cri¬ 
terion  [2,  3].  The  N  factor  is  obtained  by  one  of  a 
series  of  proprietary  or  restricted  stability  codes.  The 
code  COSAL  by  Malik  [72]  is  in  widespread  use 
in  the  U.S.^ 

We  have  performed  a  detailed  comparison  of 
amplitude-growth  curves  {N  factor  versus  arclength) 

^  COSAL  was  provided  by  NASA  Langley  Research  Center 
for  our  study. 
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obtained  with  different  options  of  COSAL,  our  own 
stability  code  LISA,  and  the  linear  PSE  code  [68].  An 
infinite  wing  with  a  64A010  airfoil  and  sweep  angles 
of  0^,  25^,  53^,  and  70®  was  considered  at  Ma  =  1.5 
and  an  altitude  of  55,  000  feet.^  The  modified  version 
WING  of  the  boundary-layer  code  by  Kaups  k  Ce- 
beci  [73]  which  is  provided  with  COSAL  was  used  to 
compute  the  basic  flow  for  all  stability  computations. 
The  flow  over  this  slender  airfoil  exhibits  weak  TS  in¬ 
stability  at  0®  sweep  and  is  otherwise  dominated  by 
cross-flow  instability. 


Figure  29:  Amplitude  growth  curves  for  cross-flow 
vortices  obtained  with  COSAL  (envelope  method), 
PSE  without  curvature  effect,  LISA,  and  PSE. 

Figure  29  compares  the  amplitude  growth  of  steady 
cross-flow  vortices  at  53®  sweep  obtained  with  the  en¬ 
velope  option  of  COSAL,  the  PSE  without  account¬ 
ing  for  curvature  effects,  LISA,  and  the  PSE.  COSAL 
does  not  account  for  curvature.  The  largest  devia¬ 
tion  between  the  results  is  caused  by  the  violation  of 
physical  constraints  by  the  envelope  method  which 
at  any  point  selects  the  spanwise  wavenumber  for 
maximum  growth  while  physical  solutions  have  fixed 
spanwise  wavenumber  [14].  Although  our  version  of 
COSAL  provides  alternative  strategies,  the  physical 
constraints  cannot  be  satisfied.  The  second  major 
change  is  the  strong  stabilizing  effect  of  surface  curva¬ 
ture.  (The  PSE  approach  does  not  suffer  from  the  in¬ 
consistency  between  in-plane  curvature  and  parallel- 
flow  approximation  [44,  43].)  The  discrepancy  be¬ 
tween  results  of  LISA  and  PSE  code  reflects  the  small 
effect  of  nonparallelism.  While  the  discrepancy  is 

^The  Euler  solutions  for  the  various  cases  were  provided  by 
G.  Klopfer,  NASA  Ames  Research  Center. 


small  in  this  case,  it  is  significant  for  other  airfoils 
and  the  growth  predicted  by  LISA  may  be  substan¬ 
tially  higher  or  lower  than  that  obtained  with  the 
PSE  code. 

The  comparison  of  the  different  results  for  N  fac¬ 
tors  enhances  the  critical  attitude  toward  tran¬ 
sition  predictions  and  raises  doubts  in  the  value  of 
the  extensive  database  of  N  factors  created  by  use  of 
COSAL  and  other  stability  codes. 

Other  studies  on  swept  wings  [41]  aim  at  comparing 
receptivity  and  disturbance  evolution  with  the  low- 
speed  experiments  of  Radetztsky  et  al.  [49,  50]  on  a 
45®  swept  wing  with  an  NLF(2)-0415  airfoil.^  The 
flow  over  this  airfoil  is  one  of  the  cases  where  LISA 
predicts  a  substantially  larger  amplitude  growth  than 
the  PSE  code. 


Figure  30;  Nonlinear  growth  of  an  unsteady  cross- 
flow  vortex  with  /  =  180  Hz,  (3  =  0.4,  and  an  initial 
amplitude  of  Aq  -  0.0002  at  the  neutral  point  in  the 
flow  over  the  NLF(2)-0415  airfoil. 

Figure  30  shows  the  amplitude  growth  of  an  un¬ 
steady  cross-flow  vortex  calculated  using  Af  =  3.  Ow¬ 
ing  to  the  large  growth  rate,  the  vortex  reaches  early 
saturation  near  30%  of  the  chord.  Steady  vortices 
of  the  same  initial  amplitude  saturate  downstream  of 
60%  chord  at  a  higher  amplitude  level.  Nevertheless, 
unsteady  vortices  create  a  stronger  mean-flow  distor¬ 
tion,  As  in  swept  Hiemenz  flow,  the  saturation  ampli¬ 
tude  is  largely  independent  of  the  initial  amplitude, 
however,  it  varies  with  the  spanwise  wavenumber  /?. 
A  detailed  comparison  with  the  experimental  data  of 
Radeztski  et  al.  is  forthcoming. 

5.2  Flow  over  a  Blunt  Cone  at  Ma  =  8 

The  stability  of  the  flow  over  a  blunt  cone  at  Ma  =  8 
has  been  experimentally  studied  by  Stetson  et  al.  [74] 
for  different  nose  bluntness.  The  cone  has  a  half  an¬ 
gle  of  7®,  a  length  of  40  inches,  and  a  base  diameter 

'’The  geometry  and  computational  results  for  the  pressure 
distribution  at  -4^  angle  of  attack  were  provided  by  W.  S. 
Saric. 
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of  9.823  inches.  The  stability  of  the  mean  flow  over 
the  blunt  cone  with  a  0.15-inch  nose  tip  is  considered 
here,  because  this  case  is  controversial  and  has  devel¬ 
oped  into  a  benchmark  case  for  application  of  basic- 
flow  and  stability  codes.  While  the  computational 
tools  provide  converging  results,  the  discrepancy  be¬ 
tween  stability  results  for  sharp  and  blunt  cones  in 
the  same  experimental  facility  and  the  disagreement 
of  the  blunt  cone  results  with  stability  calculations 
remain  unexplained. 

The  basic  flow  has  been  calculated  using  two  differ¬ 
ent  codes:  a  Thin  Layer  Navier-Stokes  (TLNS)  code 
developed  by  Esfahanian  [75,  76,  77]  and  the  AFWAL 
Parabolized  Navier-Stokes  (PNS)  code.  The  TLNS 
solution  has  been  computed  to  serve  as  an  accurate 
benchmark. 
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Figure  31:  Comparison  of  spatial  growth  rates  as  a 
function  of  arc  length  using  different  solutions  for  the 
basic  flow.  F  =  82.7  •  10"®.  □:  TLNS  solution,  200 
wall  normal  points.  <C>:  AFWAL  PNS  solution,  501 
radial  points.  A:  AFWAL  PNS  solution,  251  radial 
points. 

Figure  31  compares  the  spatial  growth  rates  for 
the  two-dimensional  second  mode  at  F  =  82.7  •  10“® 
obtained  with  different  basic  flows.  Although  the 
PNS  code  has  been  pushed  to  the  limits  with  501  ra¬ 
dial  points  (and  an  automatic  increase  in  the  number 
of  streamwise  stations),  there  remains  a  significant 
difference  from  the  growth  rates  obtained  with  the 
TLNS  solution.  This  difference  is  enhanced  by  inte¬ 
grating  the  growth  rates  into  N  factors.  Comparison 
of  basic-flow  quantities  shows  in  general  good  agree¬ 
ment  but  differences  near  the  edge  of  the  boundary 
layer.  Since  the  critical  layer  of  the  second  mode  is  lo¬ 
cated  in  this  region,  the  stability  results  are  strongly 
affected  by  the  accuracy  of  the  solution. 
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Figure  32:  Nonlinear  evolution  of  a  2D  second  Mode 
at  F"  =  82.7  •  10"®  for  an  initial  amplitude  of  0.0002. 


The  nonlinear  evolution  {N  =  4)  of  a  two- 
dimensional  second  mode  at  F  =  82.7  •  10“®  and  a 
small  initial  amplitude  of  only  2  •  10“"^  is  shown  in 
figure  32.  The  TLNS  solution  has  been  used  as  ba¬ 
sic  flow.  In  the  linear  case  (figure  31),  the  second 
mode  exhibits  growth  well  beyond  R  =  2800.  Here 
the  nonlinear  terms  cause  a  stabilization  starting  near 
R  ^  2500  and  decay  of  the  mode  for  R  >  2760.  In 
spite  of  the  moderate  maximum*>amplitude,  the  dis¬ 
turbance  causes  a  significant  rise  in  the  skin  friction 
coefficient  Cj  that  dramatically  increases  for  larger 
amplitudes. 

We  have  studied  various  second-mode  and  first¬ 
mode  interactions  to  identify  resonance  mechanisms 
and  drive  the  flow  into  transition.  Temporal  analy¬ 
sis  of  subharmonic  transition  at  Ma  =  4.5  [78,  79, 
80]  has  shown  a  strong  resonance  between  a  two- 
dimensional  second  mode  and  a  three-dimensional 
subharmonic  first  mode.  In  the  blunt-cone  flow,  spa¬ 
tial  analysis  of  this  mechanism  has  only  a  weak  effect. 
The  growth  of  the  subharmonic  mode  is  enhanced 
only  in  a  small  range  of  Reynolds  numbers,  and  as 
the  second  mode  grows,  it  suppresses  the  growth  of 
the  first  mode  below  the  linear  level  until  it  decays. 
To  clarify  this  unexpected  observation,  we  performed 
various  spatial  simulations  of  the  subharmonic  “tran¬ 
sition”  at  Ma  =  4.5  in  the  flat-plate  boundary  layer. 
The  parameters  were  chosen  to  match  those  of  the 
temporal  studies  and  later  varied  in  their  neighbor¬ 
hood.  We  found  that  the  lack  of  a  strong  subhar¬ 
monic  resonance  is  not  particular  to  the  blunt-cone 
flow.  In  the  spatial  case,  enhanced  growth  of  the 
subharmonic  is  limited  to  a  short  region,  too  short  to 


increase  the  amplitude  of  the  (assumed)  small  sub¬ 
harmonic  mode  to  a  level  where  it  would  strongly 
interact  with  the  two-dimensional  mode.  As  in  the 
temporal  case,  we  found  no  evidence  for  fundamental 
secondary  instability. 

Of  all  the  cases  studied  in  the  flow  over  the  blunt 
cone,  the  strongest  interaction  and  growth  was  ob¬ 
served  for  an  interacting  pair  of  oblique  second  modes 
with  F  =  82.7  •  10“^  and  an  initial  wave  angle  of 
40^.  With  initial  amplitudes  of  2  •  10“^  dX  R=  2000, 
breakdown  occurs  near  the  end  of  the  cone.  The  lon¬ 
gitudinal  vortex  components  exhibit  strong  growth. 

5.3  Heat  Transfer  on  Turbine  Blades 

This  last  group  of  PSE  studies  is  concerned  with  the 
heat  transfer  in  transitional  flows  and  the  prediction 
of  transition  in  the  highly  disturbed  boundary  lay¬ 
ers  on  the  strongly  curved  surfaces  of  gas  turbine 
blades  [81].  Various  experiments  are  available  that 
address  specific  aspects  such  as  the  effect  of  the  tur¬ 
bulence  level  on  the  transition  location  in  the  flow 
over  a  heated  flat  plate  [82],  the  effect  of  curvature 
on  transition  in  flows  over  heated  plates  [83,  84],  and 
the  effects  of  turbulence  on  the  heat  transfer  of  a  ro¬ 
tating  turbine  model  [85,  86]. 

We  have  chosen  this  set  of  experiments  to  explore 
the  feasibility  of  formulating  proper  input  models  for 
the  PSE  computation  that  describe  the  disturbance 
environment  and  permit  transition  predictions  with¬ 
out  empirical  N  factors.  This  input  model  was  inten¬ 
tionally  kept  as  simple  as  possible. 

Disturbances  can  be  introduced  into  the  PSE  com¬ 
putation  either  by  internal  forcing,  forcing  at  the 
boundaries,  or  through  the  initial  conditions.  The 
latter  method  has  been  chosen  in  most  DNS  studies 
such  as  the  work  of  Rai  k  Moin  [87]  and  was  adopted 
for  the  first  phase  of  our  studies.  In  flows  over  real¬ 
istic  bodies  this  method  reveals  severe  shortcomings 
when  the  computation  proceeds  through  regions  of 
strong  stability. 

The  input  model  should  introduce  a  proper  span- 
wise  scale  and  proper  frequencies  which  are  closely 
coupled  with  the  streamwise  scales.  The  spanwise 
scale  can  be  either  chosen  from  KendalPs  observa¬ 
tions  on  the  scale  of  the  Klebanoff  mode  or  from  sec¬ 
ondary  instability  theory  which  predicts  strongest  K- 
type  instability  for  spanwise  scales  smaller  than  the 
dominant  streamwise  scale  [4],  «  Aar/\/2.  This 

latter  scale  is  usually  smaller  than  the  spanwise  scale 
of  the  Klebanoff  modes. 

The  proper  frequencies,  in  particular  at  lower  tur¬ 
bulence  levels,  are  associated  with  unstable  modes 
of  the  Orr-Sommerfeld  equation  since  small  distur¬ 
bances  must  be  enhanced  by  instability  to  cause  tran¬ 
sition.  In  general,  these  unstable  modes  arc  waves  in 
a  certain  band  of  spanwise  and  streamwise  wavenum¬ 
bers. 


Since  nonlinear  interactions,  secondary  instabili¬ 
ties,  and  the  onset  of  transition  require  fluctuation 
amplitudes  to  exceed  some  threshold,  the  amplifica¬ 
tion  must  be  larger  at  lower  disturbance  levels.  Since 
the  total  amplification  between  branch  I  and  branch 
II  decreases  as  the  frequency  increases,  the  relevant 
frequencies  decrease  with  the  disturbance  level.  This 
shift  and  the  streamwise  distance  necessary  to  achieve 
sufficient  amplification  are  obvious  reasons  for  the  in¬ 
crease  of  the  transition  Reynolds  number  with  de¬ 
creasing  disturbance  level.  For  high  turbulence  lev¬ 
els,  the  weak  amplification  at  high  frequencies  may 
be  sufficient  to  initiate  transition.  At  even  higher  lev¬ 
els,  instability  is  unnecessary,  since  the  initial  distur¬ 
bances  are  large  enough  to  nonlinearly  interact  with¬ 
out  additional  amplification  and  cause  bypass  transi¬ 
tion. 

We  conclude  that  one  of  the  key  ingredients  in  the 
transition  model  is  a  pool  of  2D  waves  of  different  fre¬ 
quency  which  can  be  analyzed  simultaneously  or  sep¬ 
arately  for  their  interaction  with  the  Klebanoff  mode. 
The  upstream  shift  of  transition  as  the  turbulence 
level  increases  will  be  a  combination  of  the  upstream 
shift  at  fixed  frequency  and  the  earlier  onset  at  higher 
frequency.  The  broad  band  of  frequencies  present  in 
the  pool  of  disturbances  is  a  remarkable  difference 
from  the  traditional  vibrating  ribbon  experiments. 

The  simplest  model  consists  of  a  Klebanoff  mode  of 
spanwise  wavenumber  (3  and  a  TS  wave  of  frequency 
cj.  The  frequency  can  be  estimated  using  amplitude 
growth  curves  for  one  baseline  case  (grid  0  of  Sohn 
k  Reshotko  [82])  which  also  provides  an  estimate  for 
the  initial  amplitude.  The  amplitude  estimate  can  be 
replaced  as  soon  as  receptivity  coefficients  for  free- 
stream  turbulence  are  established.  The  amplitude 
of  the  Klebanoff  mode  is  less  critical  because  this 
mode  serves  primarily  to  establish  a  spanwise  scale, 
except  if  the  surface  is  concave  and  causes  instabil¬ 
ity  to  Gortler  vortices.  After  a  series  of  estimates  we 
chose  uj  =  0.05,  jS  =  0.6  and  amplitudes  of  0.3%  and 
0.1%  at  Rtx  —  1.25-10®  for  modes  (1,0)  and  (0,1),  re¬ 
spectively,  for  the  baseline  case  with  0.4%  turbulence 
level.  Since  receptivity  is  considered  a  largely  linear 
process,  both  amplitudes  are  linearly  scaled  for  other 
turbulence  levels. 

For  our  analysis  we  chose  the  minimum  of  the  skin 
friction  coefficient  Cj  or  the  minimum  of  the  Stanton 
number  St  as  an  operational  definition  of  the  tran¬ 
sition  location.  This  definition  necessarily  requires 
a  nonlinear  analysis  since  linear  disturbances  cannot 
cause  any  changes  in  skin  friction  or  heat  transfer  co¬ 
efficients.  The  runs  shown  were  performed  with  low 
truncation  {N  -  M  =  1)  until  shortly  downstream 
of  the  transition  location. 

The  boundary  layer  for  a  heated  flat  plate  with  an 
unhcated  leading  edge  was  used  to  study  the  effect  of 
the  turbulence  level.  Figure  33  shows  the  effect  of  tur¬ 
bulence  levels  of  1.2%  and  2.4%  in  comparison  with 
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Figure  33:  Stanton  number  vs.  Re^  for  different  dis¬ 
turbance  levels  reflecting  approximately  the  turbu¬ 
lence  levels  for  grids  0,  1,  and  2  of  Sohn  &  Reshotko. 


Figure  34:  Variation  of  the  Stanton  number  with  the 
Reynolds  number  Rex  for  convex,  flat,  and  concave 
plates  at  comparable  turbulence  levels. 


the  baseline  case,  0.4%,  The  transition  location  for 
0.4%  is  at  Rex  8,3  •  10^,  close  to  the  value  found 
by  Sohn  k  Reshotko.  The  shift  to  lower  values  as 
the  turbulence  level  increases  is  clearly  demonstrated. 
The  shift  is  not  as  large  as  observed  because  we  have 
maintained  the  same  frequency  for  these  runs.  Ac¬ 
counting  for  a  pool  of  waves  with  different  frequencies 
would  improve  the  quantitative  agreement.  Another 
important  consideration  is  the  streamwise  decay  of 
the  turbulence  level.  Rai  &  Moin  observed  a  decay 
from  9.6%  to  2.4%  over  the  length  of  their  first  com¬ 
putational  domain. 

To  analyze  the  effect  of  curvature  we  selected  three 
cases,  a  flat  plate  at  0.68%  turbulence  level,  a  con¬ 
cave  plate  with  a  radius  of  0.97  m  at  0.6%,  and  a 
convex  plate  with  a  radius  of  0.90  m  at  0.68%  from 
the  experimental  studies  [83,  84].  Although  the  tur¬ 
bulence  level  for  the  flat  plate  is  higher,  transition  in 
this  experiment  occurs  at  a  higher  Rcx  than  in  the 
baseline  case  of  Sohn  k  Reshotko.  The  effect  of  the 
unheated  leading  edge  has  been  neglected  since  the 
experimental  conditions  are  not  specified. 

Figure  34  compares  results  of  the  three  runs  which 
are  in  good  agreement  with  the  observations.  Re¬ 
markable  is  the  higher  transition  Reynolds  number  on 
the  convex  plate  which  originates  from  the  stabilizing 
effect  of  convex  curvature  on  the  longitudinal  vortex 
(see  figure  11)  and  delayed  onset  of  secondary  insta¬ 
bility.  The  significant  upstream  shift  of  the  transition 
location  is  not  caused  by  bypass  transition,  but  by 
strong  amplification  of  the  longitudinal  vortex.  The 
change  in  heat  transfer  occurs  in  a  still  laminar  but 
highly  disturbed  flow  without  participation  of  the  TS 
wave. 

5.3.1  Analysis  of  a  Gas-Turbine  Blade 

The  analysis  of  the  heat  transfer  on  the  first  stator 
blade  used  in  the  experiments  performed  at  United 


Technologies  Research  Center  (UTRC)  involved  all 
aspects  of  performing  PSE  runs  under  realistic  con¬ 
ditions. 

The  geometry  of  this  stator  blade  at  midspan  was 
obtained  in  the  form  of  134  nondimensional  coordi¬ 
nates  with  an  absolute  error  of  ±10“^.  The  geometry 
is  shown  in  figure  35. 


Figure  35:  The  airfoil  of  the  UTRC  first  stator  at 
midspan. 

The  inviscid  flow  in  the  cascade  was  computed  us¬ 
ing  the  PCPANEL  code.®  In  this  code,  the  maximum 
number  of  points  to  describe  the  blade  geometry  is 
limited  to  99.  An  attempt  to  extend  this  limit  was 
discontinued  since  the  closer  spacing  of  the  points 
caused  undesirable  oscillations  of  the  results.  Flow 
parameters  were  selected  for  the  test  run  R53PD1  [86] 
that  has  been  analytically  studied  at  UTRC.  The  tur¬ 
bulence  level  for  this  test  run  is  0.5%. 

The  pressure  coefficient  Cp  is  shown  in  figure  36. 
The  scatter  of  the  data  in  the  leading-edge  region  is 
caused  by  roundoff  errors  in  the  coordinates  which 

^This  code  was  provided  by  E.  R.  McFarleind,  NASA  Lewis 
Resecirch  Center. 
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Figure  36:  Pressure  distribution  on  the  upper  (suc¬ 
tion)  and  lower  (pressure)  surface  of  the  UTRC  first 
stator  blade. 


Figure  37:  Variation  of  the  Stanton  number  with  the 
arclength  s  along  the  suction  surface.  Comparison  of 
experiment,  laminar  flow,  and  PSE  results. 


alfect  slope  and  curvature  of  the  surface.  Otherwise, 
the  pressure  distribution  agrees  with  the  experimen¬ 
tal  data. 

The  boundary-layer  computations  were  performed 
with  a  modification  of  the  Kaups-Cebeci  code  to  ac¬ 
count  for  the  heat  flux  at  the  wall.  The  code  had 
already  been  modified  to  provide  curvature  and  met¬ 
ric  data.  The  agreement  of  the  Stanton  number  for 
laminar  flow  with  the  experiment  is  marginal,  one 
cause  may  be  different  choices  of  reference  quantities 
in  the  two  sets  of  data. 

Before  attempting  a  transition  analysis,  it  is  use¬ 
ful  to  gain  insight  into  the  stability  characteristics  by 
analyzing  the  pressure  distribution  and  the  boundary 
layer  characteristics  which  both  suggest  significant 
differences  of  the  stability  on  suction  and  pressure 
side.  The  short  region  of  adverse  pressure  gradients 
near  the  leading  edge  on  the  upper  surface  together 
with  Reynolds  numbers  in  the  range  of  Re  «  500 
causes  instability  to  TS  waves  in  this  region  followed 
by  a  stable  region  with  favorable  pressure  gradient  yet 
increasing  Re.  Strong  instability  will  appear  down¬ 
stream  of  X  0.5  once  the  stable  boundary-layer 
profile  has  adapted  to  the  adverse  pressure  gradient 
since  Re  w  750  in  this  region.  The  nondimensional 
frequencies  F  relevant  to  instability  in  this  region  are 
typically  of  the  order  F  ^  SO  •  10"^  and  higher  for  de¬ 
celerated  flows  which  translates  to  physical  frequen¬ 
cies  /  >  5  KHz.  The  disturbance  spectrum  given 
by  Bring  et  al.  [85]  provides  no  information  in  this 
frequency  range.  The  strong  fluctuation  at  the  rotor 
passing  frequency  of  190  Hz  and  its  significant  har¬ 
monics  are  far  outside  the  critical  frequency  range. 
A  frequency  of  8  KHz  was  chosen  for  the  analysis.  A 
linear  PSE  analysis  confirmed  a  short  region  of  weak 
TS  growth  near  arclength  s  =  0.2  followed  by  a  steep 
drop  of  the  amplitude  by  more  than  seven  orders  of 
magnitude  until  instability  resumed  near  s  =  0.6  in 
reaction  to  the  adverse  pressure  gradient. 


The  parameters  for  the  nonlinear  PSE  run  were 
chosen  as  before  except  for  the  frequency  uj  =  0.1573. 
The  turbulence  level  was  assumed  to  be  fixed  at 
0.5%.  The  result  is  shown  in  figure  37  together  with 
the  laminar  flow  and  the  experimental  data.  The 
run  was  continued  through  the  transition  location  at 
s  =  0.764  to  amplitudes  of  about  10%.  The  transi¬ 
tion  location  is  slightly  upstream  of  the  experimental 
result,  probably  because  of  overestimated  initial  am¬ 
plitudes.  The  distance  between  the  neutral  point  and 
the  transition  location  is  about  16%  of  the  axial  chord 
length. 

A  similar  study  of  the  pressure  side  revealed  a  sig¬ 
nificant  increase  of  the  heat  transfer  in  agreement 
with  the  experiments.  As  in  the  case  of  a  concave 
plate,  the  heat  transfer  is  enhanced  by  strong  Gortler 
vortices  without  interaction  with  the  TS  wave. 

Although  applications  of  the  PSE  approach  for 
transition  prediction  in  engineering  problems  are  just 
at  the  beginning,  the  results  so  far  show  that  tran¬ 
sition  can  be  predicted  with  minimal  empirical  feed¬ 
back  for  the  formulation  of  input  models  to  represent 
the  disturbance  environment.  This  capability  allows 
reliable  trade-off  analysis  between  different  formula¬ 
tions. 

In  addition,  the  PSE  approach  helps  to  close  the 
gaps  in  our  understanding  of  receptivity,  stability, 
and  transition  to  the  benefit  of  future  applications. 
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